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This work reports about the problem of identifying optimal pulse shapes amongst the different (just-threshold) 
ones, of a given duration, for stimulating nervous and muscular fibers and pacing the heart. 

Both electric and magnetic excitation are considered in the fibers case. Only electric excitation by external 
electrodes is considered in the myocardium case. 

The approach is a physical-mathematical one. Experimental or clinical results are mentioned when necessary to 
justify certain simplifying assumptions or to validate some of the results obtained from the theoretical work. 

In Part I the problem of identifying an optimal pulse shape amongst the different (just-threshold) cathodic ones of 
a given duration, in the framework of cathodic make excitation, is posed and solved. 

In Parts II, If and IV some restrictions of the approach followed in Part I are eliminated: different optimization 
criteria, biphasic pulses, anodic brake besides cathodic make excitation. 

In Part V the problem posed by the optimal pulse shapes in magnetic stimulation of fibers is posed and a solution 
is given to it. In the mathematical models of threshold dynamics constructed in the present research report, both 
for electric and magnetic stimulation, a mathematical tool known as the excitation functional (a concept not to be 
confused with functional excitation), is used. This excitation functional can be either postulated or derived. 

It can be derived by a procedure known as nonlinear modal analysis applied to suitable mathematical models of 
the electric or magnetic excitation of fibers or functional syncytia, introducing the concept of region of influence. 
How this is done is explained in Appendices 1 and 2. 

In Appendix I, after discussing the concept of interval of influence of the stimulating system on the stimulated 
fiber, nonlinear modal analysis is applied to construct a mathematical model of threshold dynamics, both for 
electric and magnetic excitation. 

In Appendix II, nonlinear modal analysis is applied to a mathematical model of the excitation of myocardium by 
pacing electrodes to study the just-threshold dynamics of the electrode-tissues-myocardium system. 


Unfortunately, the name "excitation functional" was introduced by the author, a long time ago, to mathematically 
model just-threshold pulses, in the volume conductor of the tissues, when the purpose of stimulation is to excite a 
certain element in that conductor of volume. 

Of course, functional excitation refers to a practical procedure intended to excite certain nervous or muscular 
fibers to restore, as much as possible, certain physiological functions. So, it should not be confused with the 
mathematical concept excitation functional. 

But the possibility of confusion is there, at least latent, so it seems advisable to change the name "excitation 
functional". A possible alternative could be "threshold functional". 


A sample of journal articles, conference papers and patents, whose content is directly related with the subject of the present 
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-Tahayori, B.; Dokos, S. (2012) Optimal stimulus current waveshape for a Hodgkin-Huxley model neuron", Engineering in Medicine and 
Biology Society, 2012 Annual International Conference of the IEEE: 4627 — 4630. 

-Tahayori, B.; Dokos, S. (2013) Optimal stimulus profiles for neuroprosthetic devices: Monophasic versus biphasic stimulation, 
Engineering in Medicine, and Biology Society 35th Annual International Conference of the IEEE: 5978 — 5981. 

-York, R., Kokones, S., Carlton, K., Greszler, A., (2013) Neuromodulation using energy-efficient waveforms , Patent No. 8 509 903 B2. 
-Intelect Medical Inc. Boston, MA (US), (2014) Neuromodulation using energy-efficient waveforms , Patent No. 8 805 515 B2. 

-Long H, Yang G, Ma K, Xiao Z, Ren X. (2017) Effect of different electrical stimulation waves on orientation and alignment of adipose 
derived mesenchymal stem cells, Chinese Journal of Reparative and Reconstructive Surgery.31(7):853-861. 

-Xinxin Li, Shirong Qiu, Menghui Li, Hong Guo, Liming Li (2017) A mathematical modeling research for effects of electric stimulus 
waveforms on the excitability of optic nerve fibers. 8th International IEEE/EMBS Conference on Neural Engineering (NER). 

-Banerjee, A., Nag, S. (2023). Energy-Efficient Electrical Stimulation Systems. In: Thakor, N.V. (eds) Handbook of Neuroengineering, pp 
825-850 Springer, Singapore. 
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PART I - A SIMPLIFIED MATHEMATICAL MODEL TO DETERMINE OPTIMAL PULSE SHAPES 
AND DURATIONS FOR CATHODICMAKE EXCITATION OF BIOLOGICAL TISSUES 
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Abstract- Our purpose is to determine an optimal pulse shape and duration of cathodic current for pacing —just at its 
threshold- a given target element in an excitable tissue. The analysis is restricted to cathodic make excitation. 

This current can be due to an electric field produced by electrodes, or an electric field induced by a variable magnetic 
field in the volume conductor formed by the biological tissues. 

The target element may be a single cell (such as nerve fiber or skeletal muscle fiber) or a whole network of electrically 
coupled cells (such as the myocardium or visceral smooth muscle). 

The response of the target element to the applied pulse is the all-or-none emergence of an action potential. 

We calculate the corresponding minimum of dissipated electric energy and then we study its variation with the 
duration of the pulse. As will be seen, the approach to this optimization problem explained in this first section of the 
report is closely related with the theory of matched filters in electric engineering. 


In this introductory section, we use the criterion of least electric energy dissipated in the tissues to select an 
optimal pulse shape amongst the different (just-threshold) cathodic ones of a given duration, in the framework of 
cathodic make excitation. We assume that the considered tissues behave as ohmic volume conductors. 

The possible influence on threshold of the compensatory anodic pulses, that in practice follows cathodic pulses 
for electrochemical reasons, is not studied. This problem will be studied in Part III. 

From a phenomenological standpoint the excitation process of a given target element in an electrode-tissues 
system may be described by an excitation functional (1) (2). ! 

For cathodic make excitation with a just-threshold pulse i(t) of duration tp, the criterion of the excitation 
functional (in the linear approximation) reduces to the following equality : 


| ¢(tp—1)-i(t)at =a, [1] 
0 


Here qu is a certain threshold charge and g(t) is an impulsive response function studied elsewhere (2) (3) (4) 
(Figure 1). 


Time 


Fig. 1: The impulsive response function of an excitation functional. 


' Bibliographical references in parts I, III, IV, V, and in Appendix 1 are cited by numbers between curved parentheses, like 
(1), (2), etc. In the list of references the articles appear in the order in which they are cited in the text. In parts II and in 
Appendix 2 references are cited by author name (s) and year of publication and given in alphabetic order in the list of 
references. The equations are in all cases numbered by square brackets like [1] ,[2], etc. 
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This function is positive between t=0 and t=t., where t, is the useful time up to the rheobase (defined using 
strength-duration curves) (5) (Figure 2). 
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Fig. 2: The useful time up to the rheobase tu. 


Furthermore, to apply [1] as a starting point, we suppose that at the beginning of the current pulse the excitable 
membranes of the target element are at their rest state. 


Applying the Cauchy-Bunyakowsky-Schwarz inequality (6): 


tp tp tp 
Je (tp—1)-i(t) at <| g°(t)dt-[P(e)at [2] 
0 0 0 

The inequality is an equality if and only if: i(t) =K g (ts — t) [3] 


Now from [1], for a just-threshold pulse the left-hand member of [2] must be equal to the constant: q. . Then, for 
Ip 2 
a just-threshold current pulse: | i°(t)dt= —— [4] 
0 | aa ( t) dt 
0 
But if the biological tissues can be represented by an ohmic load of resistance Rr, the dissipated electric energy 
tp 
per pulse, can be calculated as: R, | i (t) dt [5] 
0 


For a just-threshold pulse this energy will be a minimum if inequality [2] is an equality, that is if 7 (t) verifies 
equation [3]. 

Substituting [3] into [1], we determine the constant of proportionality K. 

After that, from [3] we obtain the optimal pulse shape, for cathodic make excitation, and for a given duration tp 


between 0 and the useful time up to the rheobase tu: 
2 


it) = —“—. g(t, -1) [6] 
Je? (dae 


From [4] and [5] it follows that the energy dissipated by an optimal pulse of duration tp is given by: 
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[7] 


To use formula [6] for the optimal pulse shape and the expression [7] for the minimum of electric energy 
dissipated in the tissues, both qu and g(t) must be known. 

In principle, they can be determined from the measurement of the corresponding strength-duration curves for each 
target element in its electrode-tissue system. 

To see how this can be done, let us consider a rectangular pulse of cathodic current, of amplitude io and duration 
tp. From the criterion given by equation [1], the pulse will be just threshold if and only if: 


=.) [8] 
[e(o) dt 
0 


In this formula, tp is between 0 and tu. Now, if io=io(te) is known (measuring the corresponding strength-duration 
curve), from equation [8] it follows that 


q, = lim (t, -iy@,)) [9] 

d |} 
OD area 10 
8(tp) A) f.] [10] 


for each tp between 0 and tu. 

In practice equation [10] shouldn’t be used to calculate g(t) form experimental data. 

The corresponding numerical problem has a very bad condition number, so that there is too much uncertainty in 
the estimation of g(t). 

To bypass this problem, we can use a closed form analytical formula with a few adjustable parameters to 
represent g(t), substitute it in formula [8] and adjust the parameters to an experimental strength-duration curve. 
The required formulae can be obtained from suitable mathematical models of threshold dynamics for fibers (4) (7) 
or networks of electrically coupled cells (2) (3) stimulated by external electrodes. 

Also, it is possible to directly introduce the linear version of the excitation functional for a cathodic pulse 
(Equation [1]) (1) (2), and then to use a descriptive mathematical model for g(t), with suitable free parameters. 


From equation [6] the optimal pulse shape of duration tp can be determined. 
2 


As g(t) is decreasing (between 0 and tm , see Figure 1) and convex (because Te g(t) >0), the optimal pulse 
t 


shape is increasing and convex. 

From [7] it follows that as tp increases from zero, the minimum dissipated energy per pulse decreases at its turn 
and reaches its minimum when tp=tu . 

The experiments show that the parameters of the strength-duration curves, and therefore qu and g(t), depend not 
only on the excitability properties of the membranes of the target element, but also depend on certain geometric 
and bioelectric parameters of the system formed by the electrodes, the tissues, and the target element within the 
tissues. 

For example, if the target element is a fiber, they depend on the radius of the fiber. 

Besides, they depend on the distance between the cathode and the target element, and of the size and shape of the 
active electrode if it is located close enough to the target element. 

Furthermore, the dissipated energy per pulse is proportional to Rr. But Rr usually depends on the size and shape 
of the cathode (we suppose monopolar stimulation here) and it is very sensitive to conductivity variations in the 
tissues located near the active electrode (such as the ones related with inflammatory reactions to implant 
electrodes). 

Therefore, the optimal pulse shape for a given duration can be modified, in the same electrode-tissues-target 
element system, due to the variation of the parameters of the system. 
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A new digital stimulator for threshold research and synthesis of optimal pulse shapes was developed, considering 
these results (8) (9). 


CONCLUSIONS 
(a) For a given duration, and a just-threshold pulse, the optimal shape is given by equation [6]. 
(b) The energy dissipated by an optimal pulse shape decreases as the pulse duration increases. 


(c) Besides being increasing and convex, there doesn’t seem to be other general properties of the optimal pulse 
shape. For each target element in its corresponding electrode-tissues system, the optimal pulse shape should 
be determined from strength-duration measurements. 


(d) In the case of implant electrodes during chronic pacing, the variations in the parameters of the strength- 
duration curves for the same electrode-tissues-target element system could be so significant that its influence 
on the optimal pulse shape can’t be neglected (for example, when the distance between the electrode and the 
target element varies). 


(e) The present approach, related with the theory of matched filters in electric engineering, is simple and 
straightforward. But it is suitable only for the criterion of minimum dissipated electric energy per pulse, in 
ohmic tissues. In other cases, the more powerful methods of functional analysis must be used (1). 
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PART II - A LESS SIMPLIFIED APPROACH TO OPTIMAL PULSE SHAPES FOR 
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Abstract-The posed problem is to find an optimal cathodic pulse shape for chronic autonomous pacing at threshold a 
given excitable tissue, in the framework of cathodic make excitation. This problem can be approached from two 
standpoints: minimizing electric energy dissipation in tissues or minimizing battery drainage. 

An excitation functional is introduced and defined as the temporal convolution of an impulse response function with a 
cathodic pulse shape. Then a given cathodic pulse shape is just threshold if and only if the value of the excitation 
functional equals a threshold charge. Both, impulse response function and threshold charge can be experimentally 
determined from measured strength-duration curves and some additional information. Electric energy dissipated in 
the tissues (ohmic resistance) and in the whole system (in series combination of ohmic resistance of electrode’s cable, 
tissues resistance and polarization resistance of electrode-tissue interface - given by an Abel’s operator) is described 
by the corresponding dissipation functional of the pulse shape. 

The variational problem of finding the cathodic pulse shape that minimizes one or another dissipation functional 
meanwhile the excitation functional remains at its threshold value (or fixed over it) is solved by methods of functional 
analysis. The obtained pulse shape is a monotonic increasing function of time, in good accord with known 
experimental facts for cardiac pacing. Bioelectric properties of tissue’s target region -among other variables- 
determine optimal pulse shapes. 

Optimal cathodic pulse shape is a property of the electrode-tissue system. For the same electrode-tissue system, it can 
vary during chronic pacing. But if wished, optimal pulse shapes can be found from measurements with rectangular 
pulses. 


Keywords: autonomous chronic pacing at threshold of an excitable tissue, optimal cathodic pulse shapes, autonomous 
chronic pacing of the heart, electric energy dissipation in tissues, battery drainage, ohmic resistance of electrode’s 
cable, tissues resistance, polarization resistance of electrode-tissue interface, dissipation functional, excitation 
functional, restricted variational problem, functional analysis methods. 


During the first stages in design and construction of implantable cardiac pacing systems, the problem of 
identifying pulse shapes and durations that minimize the charge drainage from the batteries was posed (Angelakos 
and Torres, 1964; Roy and Wehnert, 1971; Fiandra, 1985, chapter 7). 

Now, this is not a pressing problem as it was before, but minimizing the energy dissipated in the tissues in 
conditions of acute and chronic stimulation continues to be a goal to achieve. 


The problem to be solved here, from a physical-mathematical point of view, is to find an optimum cathodic pulse 
shape for chronic autonomous pacing at threshold a given target element in an excitable tissue. 

The target element may be a nerve fiber or skeletal muscle fiber (Rattay, 1990) or an electrical syncytium such as 
myocardium or electrically coupled visceral smooth muscle. The response of the target element to the applied 
pulse is the all-or- none emergence of an action potential. 


In principle the optimization problem can be approached from two standpoints: minimizing electric energy 
dissipation in tissues or minimizing battery drainage. In this research, to simplify the mathematical analysis, we 
work with a linear model of the electrode-tissue system relating the voltage between an active electrode (cathode) 
and an indifferent one (the anode), with the electric current injected in the biological tissues through the 
stimulating system. 


> This part of the report is based on R. Sudrez-Antola (1994) “Optimal pulse shape for pacing excitable tissues” (Physics in 
Medicine and Biology, Vol. 39A, Part 1, p.432) 
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From a phenomenological standpoint the excitation process of a given target element in an electrode-tissues 
system may be described by an excitation functional (Fiandra, 1985, chapter 8; Suarez-Antola, 1994). In this 
research, also to simplify the mathematical analysis, we use a linear approximation to the excitation functional. 


To solve the optimum problem, we introduce two functionals: an excitation functional w [i] and an energy 
dissipation functional gli | : 


The excitation functional: For a cathodic make excitation with a pulse i(t) of duration’, the excitation 


functional gives a necessary and sufficient condition to identify a just-threshold cathodic pulse. 


This excitation condition, in the linear approximation, reduces to the following equality, where Q, represents a 


certain threshold charge, and G(t) is an impulse response function: 
tp 
vli]=(g.i)= [Gl -1)-i()-dr=@, 1] 
0 
ty 
Here g(t) = Ge ie t) and ( gf ) = | g(t)- f (t). dt is an internal product in the space clo.r, | of continuous real 
0 
functions in the interval lo, t, | ; 


The impulse response function G(r) is positive between t =0 and t=1t, where f, is the so-called useful time up 


to the rheobase (see Figure 1) 


Time 


Figure 1: The impulsive response function of a linear approximation to the excitation functional (Suarez- 
Antola, 1994), 


To apply Equation [1] as a starting point, we suppose that at the beginning of the current pulse the excitable 
membranes of the target element are at their rest state. 

Both, the impulse response function, and the threshold charge can be experimentally determined from measured 
strength-duration curves and some additional information. 


The energy dissipation functional: The energy dissipation functional gli can be obtained from the power 
dissipated by a cathodic pulse in the active electrode interface, in the electrode leads and, in the tissues, neglecting 
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the counter-electrode effects. If V(t)is the voltage between the input end of the active electrode lead and the 


lp 
indifferent electrode, then we put: gli | = | V(t)- i (t). dt In the linear approximation (Suarez-Antola, 1994): 
0 


V()=R, -i(e)+ | K(¢—w)-i(u)-du a 


The ohmic resistance R = R, +, is the sum of an ohmic lead resistance R. and the tissues resistance R,. This 


last one, in the conditions of electric stimulation can be taken as ohmic also (Schaldach y Furman, 1975; Fiandra, 
1985; Rattay, 1990). 


The polarization resistance of electrode-tissue interface is described, in the linear approximation, by the 


t 
convolution term | K(t — u): i(u)- du. 
0 


Then, for a cathodic make excitation with a pulse i(t) of duration ¢,, the energy dissipation functional gli | may 


be written, in the linear approximation: 
oli]=R, -[?()-ar+ faor{ fat ong dit 
0 0 oO 


In inner product notation we have: 


pli]= R, -(i,i) + (Ki,i) [3] 
Here, by definition, the linear operator K verifies: 
t 
K i(t)= [x¢ —u)-i(u)- du [4] 
t m 
According to De Boer, R. and A. van Oosterom (1978) the kernel of this operator is given by K (t) = K, {<| , 
0 


where K,,7, and mare positive parameters (O < m <1). It is an Abel’s operator. The properties of the interface 


between the active electrode and the tissues can vary during chronic stimulation (Schaldach, M. y S. Furman, 
1975; Irnich, 1985; Rattay, 1990) so we can expect some variation in the above introduced parameters. 


The energy dissipation functional gli | is positive definite because both dissipation terms, RF, (i, i ) and (Ki st ) ‘ 


are strictly positive if the cathodic current pulse is not identically zero. 


The optimum problem: To find a just-threshold cathodic pulse shape, (t), continuous in the time 


interval lo, t, | , such that: 
(1D wlio | = Q,, (just-threshold condition) while 


(ID gli, | < gli, +6 i for every other pulse shape i) + Oi such that wii, +6 i =Q, (so iy + O07 is any other 
just threshold pulse shape). 


Notice that, due to the linearity of the excitation functional, from y li, | =Q,and y li, +6 i = Q, it follows that 


for every variation 67 : 
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y|oi|=0 [5] 


From Equation [3] we derive: 
lis + 5i]—oli,]=2-R,, «(ip 57) + (Kiy, 51) + (ip, K5i) + 0(6) [6] 


Here AS i ) =R,- (5 i,O1 ) + (6 i, KOi ) is always positive if Oi is not identically zero. 
So, if we find a current just threshold pulse 1, (t ) such that for every perturbation 67 that verifiesy [5 i | =0,we 


have2-R,- (ig, di) + (Kiy, di) + (ig, K6 i) =0 , then i, (t) will be an optimum pulse shape. 


Now, the continuous operator K hasan adjoint operator K* (Cotlar and Cignoli, 1974; Griffel, 1983). 
As consequence: (ig, K6é i) = (Kg, di) [7] 


Taking [7] into account we derive: 2-R,, - (i), i) + (Kiy,5i) + (io, K6i) =2-R, (Ai, 5i) 
Here, by definition, the operator A verifies: 


A=I+ 


oa (K+ k°) [8] 


p 


As usual, J is the identity operator. 


Considering that, from Equation [1] it follows y [5 i| = ( g,0l ) and considering Equation [8], the optimum 
problem can be recast as follows: 
To find 1) (t) such that: (Aig, Oi ) = 0 for every perturbation that verifies ( gol ) =0 


Results and Conclusions 


If two linear functionals like (Aig, O i) and ( g,0 i) have the same null set, the functions Ai, (t) and g(t) must 


be proportional (Cotlar and Cignoli, 1974): 
Ai,(t)= A- g(t) [9] 


Equation [9] is a convolutive Volterra’s integral equation of the second kind, so it always has a unique solution 
for every continuous function g(t) in lo, 1 | (Moissiiewitsh, 1977). 


If A-lis the inverse operator: Ie (t) A: A’ g(t) [10] 


Now, as 1, (t) must be a just-threshold pulse: 


vip |= (g.i,) = Q, [11] 
Substituting [10] in [11] we obtain: 
i, = 7S - Ag [12] 
(g,A g) 


This is the optimum threshold current pulse that we were searching. 
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It can be shown (De Boer and van Oosterom, 1978) that the effect of the polarization of the interface can be often 
neglected. In this case the shape of the optimum pulse is determined mainly by the function g(t) = Ge a t). 


As will be studied elsewhere, known experimental results (Angelakos, E., and J. Torres, 1964; Roy and Wehnert, 
1971; Klafter, R., and L. Hrieben, 1976) agree with this prediction. 


The experiments show that the parameters of the strength-duration curves, and therefore the threshold charge and 
the impulse response function depend, besides the threshold properties of the membranes of the target element, of 
certain geometric and bioelectric parameters of the system formed by the electrodes, the tissues, and the target 
element within the biological tissues. 

For example, they depend on the distance between the cathode and the target element, and of the size and shape of 
the active electrode if it is located close enough to the target element. 


The ohmic term in the dissipated energy per pulse is proportional toR, =R, + Rr. 


But A; usually depends on the size and shape of the cathode (we suppose monopolar stimulation) and it is very 


sensitive to conductivity variations in the tissues located near the active electrode, such as the ones related with 
inflammatory reactions to implant electrodes (Fiandra, 1985). 

Consequently, the optimal pulse shape for a given duration can be modified, in the same electrode-tissues-target 
element system, due to the variation of the parameters of the system. 


To end: 


(f) Optimal pulse shape is a property of the electrode-tissue-target element system. For the same electrode-tissue- 
target element system, it can vary during chronic stimulation of myocardium, nerve, or skeletal muscle. 


wa 


In principle, optimal cathodic pulse shape could be found at each moment from measurements of the strength- 
duration curves and from a characterization of the polarization resistance of the active electrode interface, 
employing rectangular pulses of controlled electric current. For a given duration, and a just-threshold pulse, 
the optimal shape is given by equation [12]. 


(g 


(h) In the case of implant electrodes during chronic pacing, the variations in the parameters of the strength- 
duration curves for the same electrode-tissues-target element system could be so significant that its influence 
on the optimal pulse shape can’t be neglected (for example, when the distance between the electrode and the 
target element varies). 


(i) A nonlinear analysis of the optimum cathodic pulse shape remains to be done. 
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PART III - OPTIMAL BIPHASIC PULSE SHAPES FOR FUNCTIONAL ELECTRIC 
STIMULATION OF BIOLOGICAL TISSUES 


Roberto Sudrez-Antola 
OMNIA Sciences/Technologies/Services 
Montevideo, Uruguay 


Abstract-An analytical approach to threshold problems in functional electric stimulation and pacing is proposed, 
framed in the concept of excitation functional. This functional can be applied to nerve, muscle, and myocardium 
stimulation by external electrodes. An optimal shape for cathodic make excitation by means of a biphasic pulse is 
found, using the criteria of minimum energy dissipated in biological tissues and total charge compensation between 
the excitatory cathodic and the compensatory anodic phases. The method can be further developed and applied to 
other threshold problems in functional electric stimulation and pacing, including biphasic anode break excitation. 


As we have seen in Parts I and II of this report, electrical stimulation of excitable biological tissues by working 
electrodes, located near the target elements in which action potentials must be produced, pose several interesting 
problems from an engineering standpoint. 

Reference (1) reviews the state of the art in relation with the design of efficacious and safe protocols for 
functional electric stimulation, both from an electrophysiological and clinical point of views, considering the 
electrochemical aspects of charge injection and the mechanism of damage to tissues and electrodes. 

Biphasic current pulses are needed. 


Usually, a cathodic half pulse of rectangular shape is followed by an anodic one (rectangular or exponentially 
decreasing), often after a short delay time. 

The charge compensation thus obtained at the electrode-electrolyte interface, allows pacing the target tissue 
without a significant damage to the biological tissues located near the working electrode or the counter electrode 
(necessary to close the circuit) (1), (2). 


The purpose of the present paper is to contribute to the design of stimulation protocols that, besides being 
efficacious and safe, may be considered as optimal from an engineering point of view. 

Minimum charge drainage and transfer, minimum energy dissipation in biological tissues or in the whole working 
electrode-tissues-counter electrode system, as well as optimal durations of the cathodic and anodic components of 
the biphasic pulses and of the delay time between the cathodic and the anodic parts are some reasonable design 
aims that may be imposed for just threshold electric current pulses (3). 

So, a threshold condition that allows us to classify biphasic pulses as below, above, or just at threshold is needed 
to discuss and solve optimization problems. After that, a safety factor must be introduced to apply the obtained 
results in clinical practice. 

Thresholds for different pulse shapes can be measured in true electrophysiological experiments in realistic 
settings. Complementary information may be produced by digital simulation using standard computing codes that 
allows us to consider the effect of an external working electrode in a target fiber or syncytium with varying 
geometries. 

But as our intention is to develop an analytical approach to optimization problems, we will use the excitation 
functional as main mathematical tool. 

In a previous work a simplified version of the excitation functional, (taking into account membrane 
accommodation processes but assuming that activation processes were instantaneously relaxed to equilibrium 
with membrane voltage), was used to study just threshold monophasic cathodic pulse shapes that minimize the 
energy dissipated in biological tissues (4). 

As a more complex and realistic example of the proposed method, we discuss here cathodic make excitation with 
biphasic pulses in which the anodic pulse achieves total charge compensation. 

We perform this task using a theoretical framework grounded in a more complete version of the functional, that 
also allows to study both one-phase and biphasic anode break excitation. 
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The idea of an excitation functional is implicit in the early classical one and two factor models of excitation (due 
to Monnier, Hill, Rashevski, and generalized later by FitzHugh, although only for uniformly polarized 
membranes) (5) (6). 

However, it seems that the first explicit proposal of the functional was done in a contribution to the 1994-World 
Congress in Medical Physics and Biomedical Engineering (3). 

Here we postulate the functional and its main properties. But it is possible to give to it a rigorous biophysical 
foundation grounded in equations of membrane dynamics and the modern tools of nonlinear science [7]-[9]. 


THE “EXCITATION FUNCTIONAL” 


Definition of the excitation functional 


In the following we consider in general the electric stimulation of nerve and muscular fibres, as well as 
myocardium and other electric syncytia. 

For each electrodes-tissues-target element system, we introduce a dimensionless impulsive response function 
G(t) that characterizes an equivalent charge transfer response to an applied electric current pulse i(r) . 


A formula for G(t) and its description will be given in the next section. 


From now on, cathodic pulses will be taken as positive functions of time, and anodic pulses as negative time 
junctions. 
If the excitable membranes of the target element are in their rest state when ¢ <Oand the pulse begins 


t 
whent = 0, let us construct the following equivalent charge: de G ) = [Ge oe u)i(u)du [1] 
0 


Now, let us consider a pulse i(r) and let us take the maximum value of 4, (1 ) forO<t. 
This procedure makes a correspondence between each function i(t ) and a real number: this defines a 


functional F [i | ; 


There is a positive threshold value q,7 such that if a pulse i(t) verifies F[i]< 4,7 it is below threshold (the action 
potential doesn’t emerge in the target tissue), if F[i]> Yer it is above threshold (an action potential emerges in the 
target tissue), and if Fi] = q,7 it is a just threshold pulse. 

While the convolution that defines g,(r) is linear operation, the excitation functional F[i] is a nonlinear 


correspondence between functions and numbers, due to the appearance of the maximum of q,(r)in its definition. 


The impulsive response functions. 


Then the impulsive response is: 


Here 4,,4,,4, are kinetic parameters of the electrode-tissues-target element system. 

(With the opposite sign convention for electric current pulses, —G(t) must be used instead of G(v) ). 

The parameter 1, is related with the activation processes in the excitable membranes, and it is at least one order of 
magnitude greater than both, 2,and7,. When the distance between the electrode and the target tissue is great 
enough in comparison with the space constant of the tissue, and the accommodation processes are fast enough, A, 


and 1, are complex conjugate numbers. If not, they are real numbers and 4, always smaller than J, . 
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Heaviside’s unit step function H(t) that appears in Equation [2] is, by definition, equal to 0 if t<Oand to 1 if 


t>0. 
Figure 1 shows a typical impulsive response function, when both parameters 2, and 1, are real numbers. 


0) 
7) 
© 
fo} 
2. 
7) 
@ oT 
[oa 
@ 
2 
2 
S 
2 
£ 


Fig. 1. Sketch of an Impulsive Response Function. 


G(t) is positive between 0 andr, , and negative after. The growing of G(r), from G(0)=0 up to G(ty ) =1, is due to 


u? 
the activation processes in the excitable membranes of the target tissue. 
The negative part after r,, is due to the accommodation processes in those membranes. 


Considering the difference in the time scales of activation processes, membrane time constant and 
accommodation processes, the following approximate formulae are derived from [2]: 


1 Pa 
ea pe [3] 
" A3 (A, +Ay) 

bite ein [4] 

(4,-4)  % 

t= 24 [5] 


If the activation process is instantaneous: 2, =0o Then formula [2] may be simplified: 
A A 
G(t)= 2 (evot)__ et Nt a(t 
(") porn rare ).a() (6] 


When the accommodation processes are slow enough, 4, is at least one order of magnitude smaller than A, , and it 


is possible to make the approximation A, =0. Then: G(t ) =e"! H (1) [7] 


It is possible to express 2, as a well-defined function of the membrane parameters of the target element (including 
time and space constants), distance between the working electrode and the target tissue, and electrode’s radius (7) 
to (9). 

Recent applications of [7] to the study of time constants for cathodic make stimulation of nerve and muscle fibers 
as well as electric syncytia like myocardium, may be found in (10) (11). 


Application of the excitation functional to a cathodic pulse of constant amplitude: strength-duration curves for 
cathodic make excitation with monophasic current pulses 


Let us consider a rectangular pulse of duration, and positive amplitude /, : i(t) =I, [H (r)- H (r —t, )] [8] 
From [2] we obtain: 
r dq..(t) 
If r<z, then: q,..(1)=1,[G(u)du [9a] So: aes 1,.G(t) [9b] 
0 
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Is positive if0 <1 <t, and it is negative ifr >r, (Fig.1). 

y dq, .\t 
If ¢, <¢ then: declt)=1, [G(u)du [10a] So: ae () 1,(G(t)-G(t-1,) [10b] 


t-t, 


Is positive if G(t)> G(t-t,), itis zero in case of equality and negative if the inequality is reversed. 


For fixed pulse amplitudes, the maximum value of 4, (t) is attained whenr, =+1,. The minimum threshold 


t, 

amplitude (rheobase) 1, 2, then verifies: Le Rb: J G(u)du = Ger [11] 
0 

When 1, <1, the maximum of q, . (t) is attained after the end of the pulse, for an instant t such that (Fig2): 


G(?)=G(F -1,) [12] 


® 
” 
c 
a 
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Fig. 2 Sketch of the graphical meaning of equation (12). 


This is due to the increasing (between 0 andt,, ) and decreasing (aftert,, ) phases of G(t) 
If ty =0 (instantaneous activation of the membrane) the increasing phase disappears and the maximum of 


(t)is attained always at the end of the pulse. 


Solving [12] with [2], and considering the different time scales to make some approximations: 


7 A, — A, \—-e*" 
potas) : z A CS , wy ts os 


The function ¢,(r,)decreases from, to 0 when f, increases from: 1, =0 to t, =t, 


The excitation functional gives the following analytical formula for pulses of durationt, and threshold 


VeT 
amplitudes 7, (strength-duration curve): loz (1, ) ~prelh) * [14] 
| Gu \du 
é(t,) 
The limit of t../..7(t,) when t. tends to Ois the so called cathodic limit threshold charge: Q..7.o 
It is experimentally determined from the measured strength-duration curves. 
From (14) we find that: Gee =Q.70 [15] 


From [11] we obtain this formula for the cathodic time constant of the system: 
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tog =—-— = | Glu)du [16] 
0 


Thresholds for rectangular biphasic pulses with delay between the cathodic and the anodic parts 


Let us consider a biphasic pulse i,(t) formed by the superposition of two rectangular component pulses. The 
first one,i,(t) cathodic, has positive amplitude /,and durationt,. The second one, i,(r)anodic, has negative 
amplitude /, and durationr,. There is a time delay +, between them. The duration of the biphasic pulse is 


t.+tg+ty =t, 


The equivalent charge is the sum of two contributions: q.(t)= Vee (1)+ Geax (r) [17] 


G(t —u)i,(u).du [18] 


G(t—u)i,(w)du [19] 


dealt)=Oif O<t<t, +tq, So that in this time interval ¢,(t)=4¢,,.(t) 


After that, the negative anodic equivalent charge ¢q,, (r)adds to the positive equivalent charge and tends to 
diminish q,(t) If ty <e,(t,), the maximum of q,(t) will be less than Ie (t. + ¢,(t,))and to just reach Jer the 


amplitude of the cathodic pulse needs to be higher, raising the threshold. 
If t, >¢,(t,)this effect disappears: the anodic compensating pulse doesn’t interfere with the threshold of the 


excitatory cathodic part. 


BIPHASIC PULSES WITH CHARGE COMPENSATION AND LEAST ENERGY DISSIPATION IN BIOLOGICAL TISSUES 


Here we are going to apply the results obtained here to an optimization problem: find the shape or shapes 
i, (t)=i.(t)+i,(t) of just threshold biphasic pulses that produce a minimum of dissipated energy in ohmic 


biological tissues. 


ty 
In quantitative terms, we seek to make the quadratic functional: fi » (t)dt [20] 
0 
a minimum maintaining the excitation functional at its threshold value: F li, | = Ger [21] 
ty fs tottg+t, 
Note that: fi,’ (Qdt=[ire(t)dr+ | iPa(t)dr [22] 
0 0 t.+ty 


Besides the just threshold condition [21], we introduce two further restrictions. 
The first one is that 1, should be greater than f,, so that we could expect that the anodic part doesn’t interfere with 


the threshold of the excitatory cathodic part, even if this last one doesn’t have a rectangular shape. 
In this case: F li, | =F li. | Then the threshold condition must be applied only to the cathodic part of the pulse. 


The second restriction is complete charge compensation: if i, (t)is the optimal just threshold cathodic shape for 


minimum energy dissipation, the anodic shape i,(t) should minimize energy dissipation while verifying: 


t, tttg tty 
Jicr(Qdt+ — fi,(Qdt =0 [23] 
0 


t.+ty 
Let us assume that dec(t), given by [18], takes the threshold value q,, for t=t,+¢, In this case, applying 
Cauchy-Buniakowski-Schwartz (CBS) inequality (12), we have: 
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[a . 
qrer = [Gl +E. - u)i, (u)du < 
. [24] 
t,+&, t, 
( f G?(w)duy(fi*e(u)du) 
& 0 
The minimum of dissipated energy is attained when the inequality becomes equality, and this happens for 
O<t<t, if and only if: i.(t)=i..7(t)=K.G(t, +E, —t) [25] 


The proportionality constant is calculated substituting [25] in the just threshold condition: 4... (1 ct & ) = er 


. VeT 
We obtain: [26] 
i) G?(u)du 
Eo 
t.+ty +t, 
Now, once i,,,(t)is determined, Q, = i (t) dt is fixed by [23]. 
t+ty 


But if 1(t)=1 for every? , then applying a second time CBS inequality: 
2 


te +tg +t t.+ty +t 
f \(t) i,(t) di} <t, J i>, (t) dt 
t.+ty t.4tg 
This can be equality, and the energy dissipated in the tissues during the anodic phase is a minimum, if and only if 
the anodic pulse is a rectangular one: La (r) =K, {H(t —(t, +tg ) = H(t —(t, +tg +t, )} [27] 
a _ Qa 
The restriction [23] gives: Ky = : [28] 


Thus, given the impulsive response function of the electrode-tissues-target element system, the durations of the 
cathodic and anodic phases, and a suitable lower bound for the interval of time between these phases, the shape of 
the just threshold biphasic pulse that minimizes the energy dissipated in the tissues was derived from the 
excitation functional and the well-known CBS inequality. 
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PART IV - GENERALIZED EXCITATION FUNCTIONALS AND OPTIMAL PULSE SHAPES 
FOR THE ELECTRIC STIMULATION OF BIOLOGICAL TISSUES 


Roberto Sudrez-Antola 
OMNIA Sciences/Technologies/Services 
Montevideo, Uruguay 


Abstract-Threshold problems in the electric stimulation of nerve, muscle, and myocardium, using external electrodes, 
may be studied from a theoretical standpoint using the excitation functional. As shown previously, the linear 
approximation to the convolution that appears in the formula for the excitation functional and the theory of linear 
matched filters can be applied to the determination of optimal pulse shapes for biphasic pulses employed in functional 
electric stimulation and cardiac pacing. Here this approach is extended using Volterra series to obtain higher order 
approximations to the excitation functional. The determination of optimal pulse shapes is reconsidered in this wider 
framework, now using the methods of the calculus of variations to pose and solve the corresponding nonlinear 
optimization problems. The pulse shapes thus obtained are compared with those obtained using the linear 
approximation to the convolution. 


Key words: Optimization theory, electrical excitation thresholds of biological tissues, mathematical modeling, 
electrical stimulation, Excitation Functional, Volterra series, functional analysis. 


Electric stimulation of excitable biological cells (muscle or nerve fibers) or functional syncytia (like myocardium) 
by a working electrode located near the target element in which an action potential must be produced, is a well- 
known threshold phenomenon. 


From the standpoint of a black-box approach, the stimulation process may be described by a correspondence 
between each applied electric current history and a binary variable A that takes the value 0 if the stimulation fails 
(no action potential is produced) and the value 1 if it succeeds. 


From the point of view of a state space approach, the field of external electric currents polarizes the membranes. 
This polarization produces typically non uniform patterns of membrane potential and ionic channel activation or 
inhibition. To elicit a propagated action potential, the state of the cell or electric syncytium membranes, driven by 
the external electric current pulse, must reach a threshold state. 

In the state space of an excitable system, it is possible to distinguish, at least approximately enough ([1], [2]), a 
decaying set, an amplifying set, and a threshold manifold. The rest state of the excitable membranes is a steady 
state point located in the decaying set. If at the end of the external pulse the state of the target cell or syncytium 
belongs to the decaying set, the electric stimulation will fail (A =0). If it belongs to the amplifying set, an action 
potential will be produced (A =1). The threshold manifold just separates failure from success. If the focus of our 
interest is the behavior of the system below and up to threshold, then only the decaying set and the threshold 
manifold need to be considered. 


For clinical applications, thresholds for different pulse shapes are measured in electrophysiological experiments 
done in realistic settings. 

Complementary information is obtained by digital simulation using computer codes for nonlinear partial 
differential equations that allows us to consider the effect of an external working electrode in a target fibre or 
syncytium with varying geometries. 

However, as in reference [3], our intention here is to develop an analytical approach to the optimization problem 
that may be posed considering the different shapes that the stimulating pulses at threshold could have, and the 
energy dissipated into biological tissues by these pulses as the figure of merit to be minimized. For this reason, we 
will use the excitation functional as main mathematical tool. 


Excitation functionals can be considered part of a black-box approach to the mathematical description of 
stimulation processes. Nevertheless, they can be derived from the nonlinear evolution equations in state space that 
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model excitation dynamics for each electrode-tissues-target element system. 


VOLTERRA SERIES AND THE GENERALIZED EXCITATION FUNCTIONAL 


Definition of the excitation functional 
The correspondence between input x(t) and output y(t)in a single input-single output (SISO) static (without 


memory) and stationary nonlinear system, regular enough, may be represented by a power series: 


y(t) = Ya, -x"(t) [1] 
n=l 


For a causal system with memory, stationary and nonlinear but regular enough, this correspondence is given by 
a Volterra series: 


co t t 
y(t) = > [dt,...fdt, -k,(t-t),....0-t, x(t, )-.. x(t, ) [2] 
n=10 0) 

The functions k,, (¢,,...,¢,,), known as Volterra’s kernels, characterize this SISO nonlinear system (4) (5). 


If the kernels are symmetrical functions of their variables, the development [2] is unique. 
When t is fixed, then y(t) is a real number and [2] is a functional: a mapping between real functions defined in [0,] 


and real numbers. 
When y(t)is considered as a function, [2] is a nonlinear operator between function spaces. 


The system is considered weakly nonlinear when it can be represented by the first two or three terms of the 
Volterra series. 


For each electrodes-tissues-target element system, we introduce a correspondence between each history of applied 
electric current i(t) and an equivalent charge transfer response q, (r). 


If the excitable membranes of the target element are in their rest state when ¢ <Qand the pulse begins 
whent = 0, let us construct the following equivalent charge: 
ao t t 
q(t) = dV fdt,...fdt, -G,(t-4,....0-1, )i(t,)-..-i(t, ) [3] 
n=10 0 
From now on, cathodic pulses will be taken as positive functions of time, and anodic pulses as negative time 
functions. Then, as was already explained in (3), we define the excitation functional both for cathodic make and 


anodic break excitation: F [i | = Maximum, {q , (t )} [4] 
Next, we introduce a positive threshold value Q, such that if a pulse i(r) verifies F[i]<@Q, it is below threshold 
(the action potential doesn’t emerge in the target tissue and the binary variable A=0), if Fli]>@Q,it is above 
threshold (an action potential emerges in the target tissue and the binary variable A =1 ), and if Fil =Q, it is a just 
threshold pulse. In principle, we could have defined the excitation functional in such a way that its threshold value 
would be the number 1, simply dividing q,(t) byQ,. 

But if we give a sound definition of the kernels in [3], Q, results to be well known limit threshold charge for 


threshold cathodic pulses of vanishing durations (3). 


t 
In reference (3) we used only the first term [Ge-1)iG)ar, which is the linear term in [3], as the simplest 
0 
approximation tog,(r)that allows the introduction of the excitation functional in order to pose and solve the 
optimization problem for biphasic pulses: first a cathodic excitatory pulse, then an anodic one for complete charge 
compensation. 
However, if our interest is in the excitation process from the rest state of cell membranes up to a threshold state, 
the system can be considered as weakly nonlinear. 
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For our present purposes we assume that it is possible to make the approximation up to the second order term: 


t 


ae lati )lh ai + 


i ot [5] 
| dt, [Gy (t—t,,1-15)-i(t,)-i(t, )-dty 
0 0 


t 
The impulsive response function for the linear approximation to q,(t), that is 4 (¢ )= Jat —u)i(u)du 1S 
0 


constructed as shown in Part II. Fig. 1. “Sketch of an Impulsive Response Function” of Part III shows the 
behavior of G,(r). Let us reconsider the meaning of the above-mentioned figure. 


The growing phase from zero up to a maximum value reflects mainly the activation of sodium channels in the 
excitable membranes. 
If we choose G,(t,, )=1, then Q; is the abovementioned limit threshold charge. 


The negative phase corresponds to the recuperation of the membrane (mainly activation and deactivation of the 
potassium channels). 
The interval of time t, , measured from the initial instant of the applied pulse, is the so-called useful time up to the 


rheobase, as shown in Part III. 
Usually, t,, is at least an order of magnitude smaller thant,. From now on, we neglect the time scale of 


activation relative to cathodic pulse durations, so we putt, ~0. Then we have: G,(ty, )~ G,(0)=1 
This being the case, for a cathodic pulse q,(t)reaches its maximum value at the end of the pulse, and not a little 


later. (This delayed maximum occurs when the finite duration of the activation processes must be considered, as 
might be necessary for very short cathodic pulses). 


Application of the excitation functional to a cathodic pulse of constant amplitude: strength-duration curves for 
cathodic make excitation with monophasic current pulses 


Let us consider a rectangular pulse of duration? and positive amplitude /, : 
i()=1,{H()-Alt-1,) [6] 

At the end of a rectangular threshold pulse of amplitude /,; : 
qelt.)= Tor fale) +r [at |Gale.t2) at 


=Q; 


[7] 


As Q,, G,(t,) and G,(t,,t,) are known if the system has been already identified, then from [7] we obtain a closed 
form analytical formula for the strength-duration curve: 
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lor = lor (t,)= 


t, : ie Ot 
fot, ba +4-0, {ia [atta 
0 0 0 
{festa 8] 


0 


tt 
2{ fas [ostintbae 
0 0 


When G,(t,,t,) can be neglected we obtain the same strength-duration formula derived in [3], as it should be: 


Or 
Lp t.)=3=—— — (9] 
[Giq)-dt, 
0 
Iflim, 9 Z.7 ‘te =Qr , and G,(0)=1, it follows from [7]: 
G (0,0) =0. Working with higher order approximations it is possible to show that if n >2, thenG,,(0,0.....0)=0. 


Thresholds for rectangular biphasic pulses with delay between the cathodic and the anodic parts 


As in reference [3], let us consider a biphasic pulse i,(t)formed by the superposition of two rectangular 
component pulses. The first one, i,(t) cathodic, has durationt,. The second one, i,(t)anodic, has durationr,. 
There is a time delay 1, between them. The duration of the biphasic pulse is: 4. +1q +1, =t,- 


The equivalent charge is the sum of two contributions: Ge (¢ ) = ec (¢ ) +dea (1 ) [10] 


qeclt)= 


G(t—t, )i.(t, dt + 


o—> 


The cathodic contribution is given by: [11] 


t t 
| dt, [Gy(t-t,,1-t))-i,(t,)-i,(t,)-dt, 
0 0 


When the anodic part does not interfere with the effect of the cathodic one, the target cell or syncytium will attain 
the threshold at the end of the cathodic pulse if at the end of the cathodic phase the cathodic contribution gq, . (r.) 


to the equivalent charge is equal toQ,. This happens when the delay is long enough in comparison with the time 
scale of membrane activation or when this last time scale may be neglected. 


BIPHASIC PULSES WITH CHARGE COMPENSATION AND LEAST ENERGY DISSIPATION IN BIOLOGICAL TISSUES 


Here we are going to apply the results obtained in Part II to the following optimization problem, already solved in 
[3] using the linear approximation to ¢,. (rt): find the shape or shapes i,(t)=i,(t)+i,(t) of just threshold biphasic 


pulses that produce a minimum of dissipated energy in ohmic biological tissues. 


If the tissues behave as an ohmic volume conductor of lumped resistance R,, the energy dissipated in the tissues 


Ze 
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ty 
is given by: Ry - fi; (t)-dt 
0 
t 


So, in quantitative terms, we seek to make the quadratic functional: Dii | - fire ).dt [12] 
0 

a minimum maintaining the excitation functional at its threshold value Fii, |= Q, [13] 

Note that: [i,2(ldr=[ie).de+  [ialt).dr Then: Dlis)= Dfi.}+ Dl] [14] 


0 0 t.+ty 


Besides the just threshold condition [13], we introduce two further restrictions. 


The first one is that tr, should be greater than t,, so that we could expect that the anodic part doesn’t interfere with 
the threshold of the excitatory cathodic part. 

In this case: Fi, |= Fli,.] 

Then the threshold condition must be applied only to the cathodic part of the pulse. Furthermore, Dii.| and Dii,| 


may be minimized independently of each other. 


The second restriction is complete charge compensation: if i. (t)is the optimal just threshold cathodic shape for 


minimum energy dissipation, the anodic shape i,(t) should minimize energy dissipation while verifying: 


t. t.ttg tty 
fir Wat + fi. (t).dt=0 [15] 
0 t.ttg 
t.+ty tty 
Now, once int (r) is determined, Q,= fi. (t).dt is fixed by [15]. 
tottg 


In this case the same argument given in [3] applies. Therefore, the anodic pulse must be a rectangular one also: 


i,(t)=K,AH(t-(, +t) -H(t-@, +ta +t,)} [16] 


Kose: [17] 


The optimum anodic phase is thus always a rectangular pulse, of duration +, that in principle could be chosen at 


will, and amplitude, given by (17), depending of the total charge transferred by the optimum cathodic phase of the 
pulse. 


Let us determine this last one. The nonlinear nature of the problem does not allow the technique applied in [3], 
but the tools of the calculus of variations can always be applied. 
Let us assume that ¢, . (t), given by (11), takes the threshold value Q, fort =t,. 


Then we must minimize the functional D[i, |= fiat, under the threshold restriction on the cathodic phase of 
0 


the pulse: 


t. 


F[i,]= Jot. —ty i, (ty Jdty + 


ie t. : [18] 
[dt, | Gy (to -ty.te -ty)-ig (ty) +i (ty )-dty = Or 
0 0 


Let us consider any cathodic pulse, referred to the optimum one i,(s): 


Pee) 
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i,(t)=i, (t)+ 2, -5,(t)+ 2 -55(t) [19] 


Here ¢,,¢5 are arbitrary real numbers and s,(t),s,(r) are arbitrary continuous functions. 
From [18] it follows that D[i,]=@(¢,,¢,)should attain a minimum value for ¢, =¢, =0 , while the excitation 


functional is restricted by the just threshold condition: ai [i ‘ | = P(e, ,€9 ) =Q, [20] 


This restriction assures that the pulse is a just threshold one. 


Applying the method of Lagrange multipliers, we obtain these two equations (6): 


= ((61,22)-4-61.62))=0 21 
&x 


In [21], k =12and¢é, =e, =0. 
After a long string of calculations, because both s,(t), s,(r)are arbitrary, we obtain: 


i,@)= Zhe, -1)+2-[a,l, tt, -wi toh} [22] 


This is a Fredholm integral equation of the second kind (7). 
If the approximation is extended to G;(t,,t,,t;) and higher order kernels, nonlinear integral equations are 


obtained. 


Solving equation [21] we obtaini,(t,4) , and substituting the solution in the threshold restriction 


Fli,.|=(e,.¢.)=@, we determine 2. So, in principle the problem is solved. 


If the integral term in [22] is neglected, we obtain the optimal pulse shape already derived in Part I, as it should 


be: 2OE ~ Gilt -1) [23] 
i ; Qr [24] 
faa 
Substituting [23] in [2]), we obtain the approximation: 
G,(t, -t)+2-— 2r 
A G?(u)-du 
in()= > Gre) [25] 


Then [25] may be substituted in [20] to find’ . 


Thus, given: (a) the first two kernels of the Volterra series to represent the equivalent charge dee(t) that 


corresponds to the cathodic phase of a stimulating pulse in an electrode-tissues-target element system, (b) the 
durations of the cathodic and anodic phases, and (c) a suitable lower bound for the interval of time between these 
phases to eliminate interference between anodic and cathodic phases, then: the shape of the just threshold 
biphasic pulse that minimizes the energy dissipated in the tissues was derived from the excitation functional. 
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The optimal anodic phase is always a rectangular pulse, independent of the shape of the optimal cathodic phase. 
This last one depends mainly of the lower order kernels (first, second, third) of the abovementioned Volterra 
series. 


The derivation of the first order kernel G, (t, )from nonlinear membrane dynamics was already given in (1), (2) and 
(8). 
The corresponding derivation of the higher order kernels G,(t,,t,),G3(t,,t).t3),..., from a state space approach 


to excitation dynamics will be shown elsewhere. 
An example of the derivation of the kernels of a Volterra series representation of a SISO such that y(t)is a scalar 


state variable described by a Riccati equation may be found in (5). 

It remains to be done a comparison of the results obtained in the present research draft with the results of sound 
digital simulations with nonlinear evolution equations and of properly designed electrical stimulation 
experiments. 
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PART V - OPTIMAL PULSE SHAPES FOR MAGNETIC STIMULATION OF FIBERS 2 


Diego Sudrez-Bagnasco and Roberto Sudrez-Antola 
OMNIA Sciences/Technologies/Services 
Montevideo, Uruguay 


Abstract-An analytical approach to threshold problems in functional magnetic stimulation of nerve and skeletal 
muscle fibers was recently proposed, framed in the concept of excitation functional. Three generations of available 
equipment for magnetic stimulation are briefly considered, stressing the corresponding pulse shape in the stimulation 
coils. Using the criterion of minimum energy dissipated in biological tissues, an optimal shape for a current pulse in 
the coil that produces a just threshold depolarization in a nerve or skeletal muscle fiber is found. The method can be 
further developed and applied to other threshold problems in functional electric stimulation. 


Magnetic stimulation of excitable biological tissues by stimulation coils, located near the target elements in which 
action potentials must be produced, pose several interesting problems from an engineering standpoint. The 
purpose of the present part of this research report is to contribute to the design of electric current pulses 
(circulating in stimulation coils) to be used in magnetic stimulation protocols that, besides being efficacious and 
safe, may be considered as optimal from an engineering point of view. 
Minimum energy dissipation in biological tissues is a reasonable design aim that may be imposed for just 
threshold electric current pulses circulating in the stimulation coil. So, a threshold condition that allows us to 
classify pulses as below, above, or just at threshold is needed to discuss and solve optimization problems. 
Thresholds for different pulse shapes in the stimulation coil can be measured in true electrophysiological 
experiments in realistic settings. Complementary information may be produced by digital simulation using 
standard computing codes that allows us to consider the effect of an external stimulation coil in a target fiber with 
varying geometries (1). But as our intention here is to develop an analytical approach to optimization problems, 
we will use the excitation functional as main mathematical tool. In a previous work a simplified version of the 
excitation functional, (considering membrane accommodation processes but assuming that activation processes 
were instantaneously relaxed to equilibrium with membrane voltage, as discussed in (2)), was used to study the 
just threshold different pulse shapes in the coil (3),(4). 

Here we postulate a simple version of the functional and its main properties for magnetic stimulation. But it is 
possible to give to it a foundation grounded in equations of membrane dynamics and the modern tools of 
nonlinear science. This has been done recently for non-myelinated and skeletal muscle fibers (3), (4). 


THE EXCITATION FUNCTIONAL FOR MAGNETIC STIMULATION 


Definition of the excitation functional 
In the following we consider in general the magnetic stimulation of nerve and muscular fibers. 
For each system composed by the stimulation coil, the tissues and the target fiber, assuming that the membrane is 


at rest prior to t=0, we convolve time derivative of the electric current in the coil “c () with a non-dimensional 
dt 


impulse response function Gm(t) to obtain a time function im(t) which has the dimensions of an electric current. 
Here, by definition, Gy(0)=1. im(t) is an auxiliary function (not the membrane current) (3), (4). 


in(=[6y 6-0) 29 ay [1] 


The derivative of the coil current appears in (1) related to the fact that the electric field gradient in each point of 


3 This part is based on the conference article by Diego Sudrez-Bagnasco and Roberto Sudrez-Antola, Optimal pulse shapes 
for magnetic stimulation of fibers: An analytical approach using the excitation functional, 32nd Annual International IEEE 
EMBS Conference, August 31 - September 4, 2010, Buenos Aires, Argentina. 
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. : dia\t : be : Coe 
the tissues is the product of HO ana a vector function of the position of the considered point in the volume 


conductor (3), (4). 
In a second step, the maximum of im(t) is taken from t=0 onwards and compared with a threshold value Itno 


(Fig. 1). 


THRESHOLD 
ia oe 
Fs sail pot ~ 


“UNDER THRESHOLD 


=. 


— 


0 t 


Fig.1. Under threshold, just threshold and above threshold time evolution of im(t) (taken from (2)) 
A history of current in the coil is just threshold when: max ‘i (¢ ) =Itho [2] 
t>0 


Irno (a characteristic parameter of the system) allows the definition of a binary variable Am, that takes the value 1 
if the stimulation succeeds and 0 if it fails (3),(4). 

This procedure makes a correspondence between each function i,(t) and a real number: this defines a functional 
Fi. | = Max ,59 fi (t)}. While the convolution that defines i, (t) is linear operation, the excitation functional F[i] 
is a nonlinear correspondence between functions and numbers, due to the appearance of the maximum of i,,(t) in 


its definition. 


The impulse response function Gr(t) 


As shown in (3), (4), neglecting the time lag in the activation of the Na* channels, the impulse response function 
for magnetic stimulation has a similar shape to the impulse response function for electric stimulation. Fig.2, 
modified from (5) shows a sketch valid for both modalities. 


1.0 


G {t) 


. * 


Time 
Fig.2. A schematic representation of a possible impulse response function taken from (4). 
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AVAILABLE EQUIPMENT FOR MAGNETIC STIMULATION 


In general, magnetic stimulation equipment involves an RLC circuit. The R and L are mainly lumped in the 
stimulation coil (which produces the varying magnetic flux in the tissues). The current in the coil for first 
generation stimulators is an exponentially damped harmonic oscillating one. Second generation stimulators 
include additional elements (such as diodes and resistors) to the original RLC topology (7). The current in the coil 
for these stimulators is an overdamped oscillation. Third generation stimulators include additional devices (such 
as IGBT) that allows to control the current duration (pulse width). Some equipment has currents in their coils such 
that an almost square pulse of induced electric field is produced in the tissues (6) (Fig.3) 


4 

AG) 

(kA) 3} 
2 


0 1 1 1 1 1 , 
-50 0 50 100 150 200 250 300 350 400 
t (us) 

E (V/m) 

proportional 

© ay 

dt 

0 


0 50 100 150 200 250 300 350 400 
t (us) 
Fig.3. Current through the stimulation coil and induced electric field for third generation equipment 
(modified from Peterchev et al. (6)). 


These almost square pulses of induced electric field allow the determination of strength-duration curves like the 
well-known threshold curves determined in stimulation by electrodes (6). 


A MODEL OF THE MAGNETIC STIMULATION SYSTEM 
The electric load seen by the stimulation equipment can be represented by a transformer whose primary circuit is 
the stimulation coil. The secondary circuit is a distributed system that can be represented by a lumped ohmic 


resistance Rr in series with a lumped electromotive force é&r (e.m.f). This e.m.f is due to the induced electric field 
caused by the variation of the induction flux produced by the stimulation coil in the tissues (Fig.4). 


Magnetic 
field B 


——— Stimulating 
Coil 


Electric w— 


volume 
field E conductor 


Fig.4. Schematic representation of the induced electric field associated with the time variation of the 
magnetic flux in the tissues due to a stimulation coil. Modified from Jarmo Ruohonen PhD thesis, Helsinki 
University of Technology, 1998. 
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The coupling between primary and secondary circuits is given by a mutual inductance M that relates the e.m.f in 
di 
the tissues with the electric current in the stimulation coil: ép =- = [3] 
é 
The electric current in the tissues ir is given in the lumped model by: lp = a [4] 
T 
tp 
The energy dissipated by a pulse of duration tp is: En = Jer (t). iy (t): dt [5] 
0 
, 1 fo 
From [4] and [5] it follows that E p= 5- | é_'(t)-dt [6] 


This is the energy dissipated in the volume conductor of the tissues by a stimulation pulse of duration tp. 

Let us consider a current through the stimulation coil that produces an induced electric field that depolarizes the 
target fiber (virtual cathode). Under these conditions, from an analysis by means of the excitation functional (4), 
in the case of a first-generation equipment, the threshold condition given by [2] should be attained for t<d (first 
depolarizing phase of the alternating current induced in the tissues). After that, as can be seen in the lower graph 
of Fig.5, the next half wave tends to cancel the effects of the first one that extends from 0 to d. 


Fig.5. Electric current in the stimulation coil for a first-generation equipment and its time derivative (the 
electric field induced in the tissues is proportional to this derivative). 


The integration time for the calculation of the dissipated energy in this case is greater that the duration of a just 
threshold pulse until the threshold is attained. 
. 2 
[as o] -dt= Diic| [7] 


From [3] and [6] we obtain: E, a 


2 te 
| 


T 0 


The parameters M and Rr are properties of the system. The electric current in the stimulation coil ic(t) may be 
changed through varying equipment settings (within certain limits, depending on the design of the equipment). 
The dissipated energy may be considered as a functional Dfi,] of ic(t). 


PULSES WITH LEAST ENERGY DISSIPATION IN BIOLOGICAL TISSUES 


As mentioned in the introduction, from a practical point of view we should try to attain threshold with minimum 
dissipation of energy in the tissues. This criterion stems from the grounded assumption that collateral effects of 
magnetic stimulation increase with the dissipated energy. 

So, the problem is to determine which current histories ic(t) are just threshold according to the excitation 
functional (that is F[ic]=Ino) and at the same time make [7] minimum. 

In Fig.6 we can see different histories of coil current that verify the just threshold condition but give different 
values to the dissipation functional Diic | ‘ 


29 


OMNIA research report Biomed. Eng. 2 Optimal pulse shapes and durations for pacing excitable 
tissues by electric and magnetic fields August 2023 


— 


Fig.6. Construction of the excitation functional for different current histories. 


In general, oscillatory wave shapes in the coil produce in a region of the target fiber (where the action potential 
could emerge) an alternate between virtual cathodes and virtual anodes. As a first requisite, the beginning of the 
coil current should produce a virtual cathode in the target region (because the other conditions being equal, 
threshold amplitudes are lower for cathodic make stimulation than for anodic break one). 

Let us consider a pulse of coil current of duration tp that produces a virtual cathode and is just threshold. 


Then: inlts)= fy (ty —v) 22). gu =, [8] 
0 


valid for OStp<tu , being tu the so-called useful time up to the rheobase, Fig.2. 


Following the same procedure as used in (8), we begin applying Cauchy-Schwarz inequality. 
2 ft t 2 
is f{ dic(t) 


This inequality becomes an equality if and only if : diclt) =K-G,(t, -t) [10] 


tp 


fou 1) He). a 


0 


For a just threshold current in the stimulation coil, left side of (9) has a fixed value Itno. The first term of right 
side of (9) is a system property and is also fixed (given certain system configuration). So, the only adjustable 
variable in this inequality is ic(t). From (7), we can diminish the dissipated energy diminishing the integral: 


[22] « [11] 


Both just threshold and minimum dissipated energy conditions are fulfilled when the inequality (9) becomes an 
equality. So [10] must be verified. Substituting [10] in just threshold condition [8] we obtain a value for K: 


I 
5 Th,0 [12] 
[Gu 4 (u)- du 
0 
Then, integrating (10), the desired optimum current pulse shape in the stimulation coil is obtained: 


ic(t)= 2 — - Gu(tp - at [13] 


tp 


[G.°(a)-du : 


Let us suppose that the stimulation coil is modeled as a pure inductance L. Then the voltage wave shape that must 
be applied by the equipment to the coil to obtain a just threshold pulse that minimizes the energy dissipated in the 
tissues is given by: 
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Vow (t)=L ciel) =L - Imo g(t, -t) [14] 
| Gy : (u ) -du 

The best way to obtain meaningful strength-duration Ane both for electric and magnetic stimulation is by 
means of rectangular pulses of electric field in the tissues (7). 
In the case of magnetic stimulation of fibers this is possible with equipment like the one described in (6). In Fig.7 
taken from (6) five different strength-duration curves for cTMS corresponding to five rhesus monkeys are shown. 
100 
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Fig.7. Strength-duration curves for magnetic stimulation with wave shapes as shown in Fig.3. 


In this case we have only five experimental values for each strength-duration curve. Under these conditions the 
identification of Gm(t) from experimental data as described in (3) should be done with a simpler model for Gm(t) 
(neglecting accommodation processes) than the one from Fig.2. 


The simplest possible model is given by an exponential function with only one parameter ts: 


cul] 


t 


7 -H(t) [15] 


Here H(t) is the unit step function. Using [15] in [13] we obtain: 


ie(t)=| 22 io } [u()-H(e-1,) [16] 


A formula like [16] may be obtained from a coarse electrical model of excitation given by an RC circuit for the 
membrane, an ad-hoc threshold condition for voltage, and an ohmic resistance for the tissues. 
However, its meaning is completely different. 


A sketch of [16] is shown in Fig.8 
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i{t) 


0 t, t 


Fig.8. Qualitative current shape in the stimulation coil for a just-threshold and minimum dissipated energy 
condition for a simplified Gwi(t). 


CONCLUSIONS 


For a given stimulation coil-tissues-target fiber system, there is an optimum current pulse in the coil that 
minimizes the dissipated energy in the tissues. In the framework of the present model taken for the excitation 
functional, this optimum current is given by [13]. 

Gm can be obtained from experimental data as shown in (3). 

In the model for Gy illustrated in Fig.2, the time lag in the activation of the Na* channels were neglected. 

This restriction is compatible with pulse durations not too short. 

For pulses short enough the threshold is not attained at the end of the initial virtual cathodic phase but at some 
time after the end of this phase. 

Even a simpler model is given by [15]. For a more realistic Gm, like in Fig.2, i-(t) is different from [16]. 


In the case of the damped current oscillation of a first-generation equipment, if threshold is not attained before the 
end the first virtual cathodic phase (first quarter of the first current oscillation), it will not be attained at all. 
Moreover, the following oscillations contribute only to energy dissipation and not to threshold attaining. So, to 
attain threshold while dissipated energy is minimized, requires special designed equipment. 

As a next step, the improvement in the reduction of the dissipated energy in the tissues for the optimum pulse 
shapes (derived in the present paper) versus standard pulses obtainable with commercially available equipment 
should be studied. This could afford useful information to assess the feasibility of developing special equipment 
for generating optimum pulse shapes in case the diminishing of collateral biological effects justifies it. 
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APPENDIX 1 


A Nonlinear Modal Analysis of Threshold Dynamics for the Electric and Magnetic 
Excitation of Nerve and Muscle Fibers 


Roberto Suarez-Antola and Diego Suarez-Bagnasco 
OMNIA Sciences/Technologies/Services 
Montevideo, Uruguay 


Abstract—Threshold problems in electric stimulation of nerve and muscle fibers have been studied from a 
theoretical standpoint using nonlinear modal analysis. A consequence of this approach to threshold 
dynamics, is the so-called excitation functional. Here the excitation functional is extended to magnetic 
stimulation of excitable nerve and muscle fibers. A unified derivation of the functional is done, for (non- 
myelinated) nerve and muscle fibers, by means of the nonlinear cable equation with a Fitzhugh-Nagumo 
membrane model and a generalized Rattay’s activating function. The identification problem of the 
excitation functional for magnetic stimulation, from strength-duration experimental data, is briefly 
considered. 


The goal of electric and magnetic stimulation of excitable cells is to produce (or to block) action 
potentials in suitable locations. From the standpoint of a black box approach, the stimulation process 
may be described by a correspondence between each applied electric current history and a binary 
variable A that takes the value 0 if stimulation fails and 1 if it succeeds (Fig.9). 


Under) J] 
threshold 0 Binary 


Above I | | 4 
threshold 
pulses : 
\4 
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Fig.9. Black box approach to electric and magnetic stimulation. 


From the standpoint of the stimulation equipment, the black box comprises the electrodes (and their 
leads) or the magnetic coils, the volume conductor of the tissues and the target elements (nerve or 
muscle fibers, etc.). So, this black box may be considered as an electric load seen by the stimulating 
equipment. The binary variable may be obtained through an electric measurement (detection of action 
potential by recording electrodes) or by external manifestations (like muscle twitches, function 
inhibitions, etc). 

For the electrical stimulation, there is already a theoretical tool, the excitation functional, introduced 
by R. Suarez-Antola (1), (2). 

It allows both a description and a prediction of the output given by the binary variable in the black 
box approach. 
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As each system composed by the electrodes, the tissues and the target elements are unique (due to 
different spatial, temporal, electrical and in general, physiological properties), the black box must be 
duly identified from suitable experimental data. 

The experimental strength-duration curves for a given system can be used to obtain the excitation 
functional for the electrical stimulation of the just mentioned system. 

The excitation functional opens a third way between a pure experimental approach and a pure 
computational approach (working with nonlinear cable equation for electrical stimulation). As an 
analytical tool (which parameters can be adjusted from real experimental data) it allows a 
mathematical formulation of threshold related problems of interest for biomedical engineering. For 
example, to find optimal pulse shapes given an optimization criterion (2). 

So, it could be of certain interest to try to extend the excitation functional to magnetic stimulation. 


The purpose of this part of the report is threefold: (a) to extend the excitation functional to magnetic 
stimulation produced by external coils, (b) to present a unified derivation of the functional both for 
electric and magnetic stimulation, and (c) to briefly discuss the identification problem of the 
functional from experimental data (from strength-duration curves). 

A unified derivation is possible because (from the nonlinear cable equation with a generalized 
activating function (4) (5)) it is the external electric field parallel to the fibers the responsible of both 
magnetic and electric stimulation. One of the consequences of this extension is to have a tool that 
allows to characterize the system and to predict the outcome of the binary variable given a certain 
current history in the stimulating coil. 

This tool can be used also to pose and solve certain engineering problems related to the system, like 
how to shape an input current pulse in the stimulating coil that is both threshold and optimum from 
the standpoint of minimizing the energy dissipated per pulse in the tissues (6). 

For the purposes of the present work, a suitable background in electric and magnetic stimulation may 
be found in (7). Further information about electric stimulation can be found in (8) to (10), and in (11) 
to (16) for magnetic stimulation. 


EXTENSION OF THE EXCITATION FUNCTIONAL TO THE MAGNETIC CASE 


Formulation for the electrode case 


The simplest formulation of the excitation functional corresponds to a single active electrode or a 
bipolar electrode and may be applied both to cathodic make excitation and to anodic break excitation 
(2). The excitation functional when the excitable membrane is at rest prior to t=O can be written in two 
steps. First, we obtain a function qg(t), making the convolution of a history of injected current ig(t) 
with a non-dimensional impulse response function Gg(t) characteristic of the system (black box). 


aelt)= [Gele—u)iplu)-a i" 


In this case, by definition Gg(0)=1. The injected current appears in (1) related to the fact that the 
electric field in each point of the tissues is the product of iz(t) and a vector function of the position of 
the considered point in the volume conductor (4). 

In a second step, the maximum of qe(t) is taken from t=O onwards and compared with a threshold 
charge Qtno. This threshold charge is not the charge injected by the electrode neither the charge that 
crosses the excitable membrane of the target fiber. 

A history of applied current is just threshold if and only if: 


max {gel(t)}=Omo (21 


t>=0 
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Fig.10 shows three possible outcomes of the convolution. 
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Fig.10. Under threshold, just threshold and above threshold time evolution of qr(t). 


Qtno , Which is a characteristic parameter of the system, allows the definition of the digital variable Ag 

that was presented in the introduction. As will be seen in section D, in this case Qrno coincides with 

the well-known limit threshold charge that can be obtained from strength-duration curves. 
Formulation for the coil case 


In the magnetic case, assuming that the membrane is at rest prior to t=0, we propose to convolve time 


derivative of the electric current in the coil cl) with a non-dimensional impulse response function 
dt 


Gm(t) to obtain the time function im(t). Here, by definition, Gy(0)=1. 


in = [y(-u) Se au, 3] 


The derivative of the coil current appears in [3] related to the fact that the electric field in a given 


point of the tissues is the product of dic) and a vector function of the position of the considered point 
dt 


in the volume conductor (4), (5), (6), (16). 
In a second step, the maximum of im(t) is taken from t=O onwards and compared with a threshold 
current Itno (Fig. 3). 
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t 
Fig.11. Under threshold, just threshold and above threshold time evolution of im(t). 


A history of current in the coil is just threshold when: 


max tiy (t)}=In,0 [4] 


120 
Ith (a characteristic parameter of the system) allows again the definition of a digital variable, in this 
case Am. It can be obtained from experimental strength-duration curves, as shown in section IV. 


OMNIA research report Biomed. Eng. 2 Optimal pulse shapes and durations for 
pacing excitable tissues by electric and magnetic fields August 2023 


A UNIFIED DERIVATION OF THE EXCITATION FUNCTIONAL 


The construction of the impulse response function and the threshold charge for electric stimulation, 
and the impulse response function and the threshold current for the magnetic stimulation of a (non- 
myelinated) fiber will follow a common procedure. This allows an easier grasp of the similarities and 
differences between the two situations. 

Let us begin with the nonlinear cable equation with Rattay’s activating function generalized to 
consider both electric and magnetic stimulation (4), (5), (6), (16). To simplify the derivation, we use 
the well-known Fitzhugh-Nagumo model for the unit membrane (17). If the (possibly bended) target 
fiber is represented by a curve of directed arc length s measured from a suitable fixed point of the 
fiber, v(r,s)is membrane voltage field and w(t,s)is the recovery variable field, both relative to their 


rest values, we obtain the following set of equations [6]: 


Ov O°v OE, (t.s 
tm = Ev tb-v?—e-v aw} as As? - 2. A ) [Sa] 
Tw ae =(v-y) [5b] 
Ot 


The time constant of the membrane is7,,, and7,,is the time constant of the recovery variable (under 
physiological conditions at least an order of magnitude greater than T,, ). 


The fiber’s space constant is 4,, . 
The parameters of the unit membrane model b,c, a, 7 are positive. 


Rattay’s generalized activating function is given by (4), (5), (6), (12), (16): 


aE, (ts) _ [4m Fe()-ie(t) electric 
agp =), 2 | 7 
Os A. Fy (s)- qiclt) — magnetic 


The functions F;,,(s) and F,,(s)give the spatial distribution of the perturbation produced on the 
excitable membrane by the electric field due to the electrodes or induced by the time varying 
magnetic flux due to the coil, and may be called geometric form factors for electric and magnetic 
stimulation, respectively. Now, the stimulation is always a localized phenomenon. Therefore, outside 
a certain interval of influence along the fiber, of length £ , the form factors may be neglected (6), (18), 
(19), (20). If we assume that until the fiber reaches a threshold state, the fields of membrane voltage 
and recovery variable take their rest values at the end points of this interval of influence, it is possible 
to make a nonlinear modal analysis of the cable equations. The value of 4 depends on the form 
factor. Ifsis measured from the mid-point of the interval of influence, we may use the Fourier 
development: 


v(t,s) =A, (t)- 1 


w(t,s)= B,(t)- | 


Fy (s)= Fy ‘y 
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Substituting [7] in the nonlinear cable equation and eliminating the spatial dependence, we obtain a 
system of first order nonlinear ordinary differential equations in the unknown mode 
amplitudes A; (t)and B, (1). Solving this system with suitable initial conditions, we obtain the mode 


amplitudes. The procedure is the same as already developed for stimulation by electrodes (18), (19). 
The only difference between electric and magnetic stimulation in this approach is the forcing term. 
So, the dynamics of the unforced fiber, after the end of the stimulating pulse, is the same if the length 
of the interval of influence is the same. 

Digital simulation of both cathodic make and anodic break electrical stimulation using up to 
seventeen mode amplitudes suggest that the first mode is the most relevant in determining threshold 
behavior of the fiber (20). 

Truncation to the first mode, uncoupled, gives in both cases, a set of nonlinear ordinary differential 
equations in the mode amplitudes corresponding to membrane voltage and the recovery variable (19). 
The study of threshold dynamics done with the abovementioned equations, with parameter’s values 
within the physiological ranges, shows that a threshold barrier may be defined in phase space 


(A, , B,) such that when the phase point reaches the barrier, an action potential emerges (6), (19). 


The under-threshold behavior can be described with enough accuracy, up to the threshold barrier, by 
the linear system: 


dA Aan ar 
8 at =a em go + Ay Pina wel) 
dB, 
a ( 1 Y ) 


This approach is the equivalent, for a fiber with a non-uniformly polarized membrane, to the classical 
discussion of Fitz-Hugh for a unit membrane (19). Fig 12 shows the simplified dynamics in the 


(A,B, ) state space. The results of the digital simulation suggest the introduction of a decaying and 
an amplifying set, separated by a threshold curve (19). 


y; 
#*~ Threshold barrier 


re 


+ 
Fig 12. Sketch of the dynamics of the system, simplified by the threshold barrier. R is the rest 
state. The decaying set is shown as a shaded area bounded by the double arrowed threshold 
curve. 


The threshold curve is composed by half-straight line taken from the threshold barrier and an orbit in 
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the decaying set that is just tangent to the threshold barrier. The construction of the threshold barrier 
can be seen in [17]. 


The linear dynamic system (8) can be recast in a matrix form: 


[9] 


F ~ 
x= A= ey -|3) [10] 
B u = 


2 


In the matrix A, the spatial properties appear in the element = los An through the non- 


2 


; . A 
dimensional number 7s . 


The other elements of the matrix A depend only on unit membrane properties through the parameters 
a, Y 2 Tm ) Ty . 


The solution, beginning from the rest state, and until the first arrival to the threshold barrier, is: 


- rd. = 
X(O)= dy’ “Fa fapie(u)- (ES, )- du [11] 


0 


From the results of digital simulation (summarized in Fig 12) the condition that characterizes a just 
threshold magnetic stimulation is (Fig.5): 


max vi *(0)}= Pp [12] 
te|0,+00) 


I Threshold barrier 


Fig.13. Graphical interpretation of Eq.(12). The unit vector nis normal to the threshold barrier. 


There are analytical formulae for the unit vector n and the distance p as functions of system’s 
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parameters (19), (20). 
We put x(t) from [11] into [12] and define: 


(a) The impulse response function: Gy (r ) =F [13] 


(b) The threshold current Tino [14] 


Then we obtain the just threshold condition for magnetic stimulation in terms of the excitation 
functional (equivalent to (3) and (4) combined): 


f d.ilu 
max [Gu(t-u) < ( du =Inno- [15] 
te[0,+20) 0 du 
The same procedure applied to electric stimulation gives: 
-T tA- 
Dp nee, 
Ono = G;\t)=— 
vi ve ‘Fry n'é, al ) ni é 


So, if the parameters of the unit membrane are the same, the only difference in the impulse response 
function (in the framework of this mathematical model) is due to different values of @ related with 
differences in the form factor. 

The activation of peripheral nerves can be studied under the following modeling conditions: the 
medium can be considered as homogeneous; the fiber can be considered as straight and unbounded, 
but the external electric field varies in a region along the fiber (12). 

In the framework of the present mathematical model, the degree of spatial localization of the 
perturbation of the peripheral nerve membrane due to the external field is given by the length of the 


interval of influence e, 
It is possible to show that an approximate formula for the time constant of the cathodal strength- 


duration curve for a nerve fiber is given by (21): t. = a [16] 
t +77 ae 

The chronaxie is proportional to ts. For a given peripheral nerve, in the electric stimulation / is 

smaller than in the magnetic stimulation case. From (16) it follows that chronaxies for electric 

stimulation should be smaller than chronaxies for magnetic stimulation. This explains the 

experimental findings reported in (13). 

Equation [16] is obtained from the simplest model for the impulse response function. A more realistic 

model of Gm (given by [13]) is sketched in Fig.14. 


10 


G ft) 


t t 
u ove Time 


Fig.14. Schematic representation of an impulse response function taken from (19). 
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It allows the calculation of the chronaxies both for cathodic and anodic stimulations. However, 
despite the analytical formulae are different, the same theoretical predictions about the behavior of 
chronaxies for electric and magnetic stimulation of peripheral nerves are derived from this more 
complete model. 


THE IDENTIFICATION OF THE IMPULSE RESPONSE FUNCTIONS 


For electric stimulation, the identification problem of the excitation functional is already studied in 
(2), (3), (6). 

For magnetic stimulation, let us consider a linear ramp of electric current in the coil. From the formula 
of the excitation functional [15] it follows that for a ramp of duration tp, that produces a suitable 


depolarization of the fiber membrane in the interval of influence, the threshold slope <n (tp) is 
t 


given by: =o (tp) =- ‘m0 [17] 
[Gy (u)- du 
0 


From (15) when the ramp duration tends to zero, we derive: 


: d., 
I7,9 = lim| t, Gv iemlte) [18] 
P 


Once Itno is known, from [16] and [17] it follows: 


Gy (,)=—2_| "mo — [19] 
Otp | dic ry (1) 
age 


This means that the impulse response functions may be obtained from experimental strength-duration 
curves for magnetic stimulation. These experimental curves can be done using equipment like the one 
described in (22). However, the derivative makes it an ill posed problem. A way out of this difficulty 
is to use analytical formulae [13] and [14] derived in this paper and adjust the parameters of [15] to 
the experimental data. 


FINAL COMMENTS 


The excitation functional is a systemic property. From the standpoint of attaining the threshold in a 
given target fiber in an external electric field, an action potential can emerge under four modeling 
circumstances: (a) the external field is uniform and the fiber crosses a region of fast variation of the 
electrical conductivity of the medium, (b) the external field is uniform and the medium can be 
considered as homogeneous but the fiber bends, (c) the external field is uniform and the medium can 
be considered as homogeneous but the fiber originates or terminates (short circuit or open circuit 
conditions), (d) the medium can be considered as homogeneous, the fiber can be considered as 
straight and unbounded but the external electric field varies in a region along the fiber. The model of 
the present paper applies to circumstances (b) and (d), because what matters is the spatial variation of 
the projection of the electric field tangential to the fiber and its time variation. This time variation is 
proportional to time variation of the current in the electrode case, and the time derivative of the 
current circulating in the working coils. An extension to cases (a) and (c) could be done. 

In the magnetic case the determination of the form factor is a difficult problem that deserves further 
study. 

However, once determined Gy(t) and Ito, the binary output response of a given target fiber to 
different pulse shapes in the coil could be predicted and contrasted with experiments, assuming that 
the system composed by the coil, volume conductor and target fiber remains unchanged. 
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An extension of the present derivation of the excitation functional for magnetic excitation, to consider 
lag effects in the activation of the excitation channels in fiber’s membrane, and myelinated fibers, both 
neglected in this paper, remains to be done. 
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APPENDIX 2 


A Nonlinear Modal Analysis of Threshold Dynamics for Cardiac Tissue Excited by 
External Electrodes 


Roberto Sudrez-Antola 
OMNIA Sciences/Technologies/Services 
Montevideo, Uruguay 


Abstract- Thresholds in cardiac pacing continue to be an active area of research, both from an 
experimental and a theoretical standpoint. Most theoretical research hinges on ab-initio numerical 
studies framed in bidomain models of cardiac tissue. Analytical results are scarce and almost always 
based on oversimplified models of the electrode-myocardium (E-M) system. 

Nonlinear modal analysis is applied to a realistic mathematical model of the excitation of myocardium by 
pacing electrodes to study the just-threshold dynamics of the E-M system. 

A region of influence of the external electrode on the myocardium is introduced to pose a mixed 
homogeneous Dirichlet- nonhomogeneous Neumann problem for the nonlinear bidomain field equations, 
describing the polarizing effect of the electrode on the excitable tissue and considering anisotropy effects. 
Then, a nonlinear modal analysis is used to transform the field equations in a system of nonlinear 
ordinary differential equations in mode amplitudes. The threshold value of a perturbation that causes 
instability is estimated by a modified version of Eckhaus method. 

A whole family of threshold space patterns of membrane polarization is identified in the state space 
constructed from the modal amplitudes. The sizes and shapes of liminal regions (that must be depolarized 
over threshold of uniform membrane excitation to start a propagated action potential) are obtained. 
Analytical expressions for cardiac strength-duration (S-D) curves are derived, both for controlled current 
and controlled voltage stimulating pulses, including closed form formulae for time constants (or 
chronaxies) and rheobases in terms of geometric and electrochemical electrode parameters, electrode- 
myocardium distance, electric properties of the volume conductor, and local distribution of electrotonic 
and excitability parameters of myocardium. 

Several experimental results are thus explained (liminal regions and the behaviour of the parameters of 
S-D curves) and some new facts are predicted (related with liminal regions, the family of E-M threshold 
states and the parameters of the S-D curves). This theoretical approach can be extended to other excitable 
tissues, such as electrically coupled smooth muscles. Besides, it can be used as a guide to the design of 
pacing electrodes. 


Most theoretical research related with thresholds in cardiac pacing, hinges on ab-initio numerical 
studies. These studies are framed in bidomain mathematical models of cardiac tissue (Janks and Roth, 
2009). Analytical results as such are scarce and almost always based on oversimplified models of the 
electrode-myocardium system. 


The physical foundations and the mathematical formulation of bidomain models in electro cardiology 
can be found summarized in Bradley Roth’s article in Scholarpedia and in more detail in the book 
“Cardiac bioelectric therapy” edited by Igor Efimov, Mark Kroll and Patrick Tchou that appears in the 
references of the present paper. 


The parameters of measured strength-duration curves (time constant, rheobase, limit threshold charge) 
show a dependence of electrode size, shape, and material (Laplantif, 1983; Irmmich, 1983), distance 
from electrode to myocardium (Irnich, 1975 and 1980), and tissues properties (Furman, Hayes, and 
Holmes, 1993). A full understanding of the interrelation between these variables is of interest, both 
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from the engineering (electrode’s design) and from the medical (practice of cardiac pacing) points of 
view. 


To generate an action potential, a suitable polarization pattern (relative to the rest state) must be 
attained in the excitable membranes of myocardial cells. In cardiac pacing, the perturbation from the 
rest state is produced by the active electrode in the excitable membranes of a near-by region of 
muscle. The rest state is stable to small enough perturbations. As consequence, an analytical approach 
to excitability must allow a derivation of threshold amplitudes of perturbation beyond which the rest 
state is unstable. 


To cope with stability problems in the large, posed by finite amplitude instabilities in fluid mechanics, 
an asymptotic method of nonlinear modal analysis was developed (Eckhaus, 1965). During the sixties 
and the seventies, it was applied to the study of stability problems in other distributed parameter 
systems, including mass transport processes and chemical reactions in continuous flow chemical 
reactors and other systems of interest from an engineering standpoint (Denn, 1975). In the nineties it 
was applied to study the just threshold dynamics of excitable tissues stimulated by external electrodes 
(Suarez-Antola, 1994; Suarez-Antola and Sicardi-Schifino, 1996; Korenko, 1997; Suarez-Antola, 
1999). 


The purpose of the present appendix to this report work is to apply the above-mentioned method of 
nonlinear modal analysis to a mathematical bidomain model of the excitation of myocardium by 
pacing electrodes. The focus will be the study of the just-threshold dynamics of the electrode- 
myocardium system. Analytical formulae for the parameters of strength-duration curves are obtained, 
as functions of the geometric and electrophysiological variables, considering the unequal anisotropy 
of cardiac muscle. 


As things stand today, it seems appropriate to make a few comments on the interactions between 
analytical models and numerical models. 

Many mathematical models in physiology are scarce in detail (from the physiologist point of view) 
but not trivial if you want to approach them from an analytical point of view. Consequently, research 
in mathematical physiology is mainly done with digital simulations. This poses a challenge: the 
efficient exploration of high-dimensional parameter spaces. 

So, a trade-off must be done between the inclusion of numerous physiological and biophysical details 
(for example, a hundred of different types of membrane channels) whose parameters are often 
measured with significant uncertainties that are sometimes transferred amplified to the simulation 
model, and oversimplified models that perhaps neglect important features (for example, the dynamics 
of adaptation variables in excitable membranes is not included together with activation and recovery 
variables). 

Standard computer packages allow to obtain numerical solutions for selected values of the parameters. 
It is possible to obtain suggestions about analytical (algebraic) results from the numerical information 
using software to plot multidimensional data and examining various projections. But if understanding 
is the aim, it seems advisable to try as hard as possible to obtain analytical results. This is the focus of 
the present work. 


Restricted nonlinear modal analysis of distributed parameters systems: an outline. 


In a generalized sense modal analysis is a procedure that allows: 
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(a) The expansion of the solutions (fields) of partial differential equations in time and space 
variables as a series of given space functions weighted by time functions known as modal 
amplitudes. 


(b) The derivation of evolution equations for the mode amplitudes (a set of ordinary differential 
equations). 


As the space functions are known, by means of modal analysis a complex space- time field dynamics 
is reduced to the study of the evolution of a representative point in the space of mode amplitudes. If 
the original partial differential equation is nonlinear, the ordinary equations for mode amplitudes are 
nonlinear also: in this case we have nonlinear modal analysis. An often-used procedure to obtain the 
equations for mode amplitudes is the application of Galerkin’s method. In general, in their evolution 
equations modal amplitudes appear coupled, even in the linear approximation. 


In a restricted sense, modal analysis is a procedure applied to the different fields that give the space- 
time evolution of a distributed parameters system with a bounded space domain B, expanding the 
fields in series of eigenfunctions of a certain suitably chosen linear operator L. The domain of this 
operator, which is a set of space dependent functions defined in B, is so defined to include the 
boundary conditions of the problem. These boundary conditions must be time independent, although 
they can be non-homogeneous. 


The abovementioned eigenfunctions (spatial modes) depend only of the space coordinates and, as 
already said, are defined in the bounded domain B. In each expansion corresponding to a given field, 
these spatial modes are weighted by unknown time dependent modal amplitudes. The modal 
amplitudes may be determined so that the series expansions represent a solution of the dynamic field 
equations verifying the initial and boundary conditions in the bounded space domain. 


The original nonlinear partial differential equations, with time and space as independent variables, are 
substituted by an equivalent system of nonlinear ordinary differential equations for the unknown 
mode amplitudes, with time as the only independent variable. The advantage of this restricted kind of 
modal analysis, compared with the generalized one, is that in the obtained evolution equations the 
modal amplitudes always appear decoupled to the linear order. 


From the initial conditions verified by the fields, a set of initial conditions for the mode amplitudes 
are obtained. Thus, as restricted modal analysis is part of generalized modal analysis, also in this case 
a complex space- time field dynamics is reduced to the study of the evolution of a representative point 
in the space of mode amplitudes. 


Often, it is possible to work with a relatively small number of mode amplitudes, after reducing the 
original high order dynamics to a low order one: this is the subject of inertial manifold theory. A 
necessary condition to be able to work, at least asymptotically, with a small number of time dependent 
mode amplitudes is that the eigenvalues of the considered linear operator be widely spaced. 


Now, there is a close connection between the classical Liapunov’s First Method to study stability 
problems and restricted nonlinear modal analysis. The First Method of Liapunov, also known as local 
linearization, or stability in the presence of infinitesimal perturbations, allows us to study strictly only 
the stability of the linear version of the dynamic system. However, with the aid of the well-known 
theorem due to Hartman and Grobman, its scope can be extended to study the stability of non-linear 
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systems in a small enough neighborhood of the steady state. But it can't cope with stability 
problems in the large (that is, the response to non-infinitesimal perturbations). However, the design, 
operation and control of cardiac pacing systems often pose stability problems in the large. The 
analytical methods of restricted non-linear modal analysis allow us to tackle some of these problems 
from an analytical point of view, as will be shown in this paper and as shown elsewhere (Suarez- 
Antola, 1994 and 1999). So, restricted nonlinear modal analysis can be considered an extension of 
Liapunov's First Method to cope directly with the non-linear terms in the dynamics. These methods 
complement both, numerical methods ab-initio, and the results obtained after applying some 
summarizing function criterion to determine a region of stability, like Liapunov's functions 
(Liapunov’s Second or Direct Method) and the like, that are often applied in theoretical discussions. 


Let us consider the following nonlinear evolution equation for the scalar field u(t a); as a simple 


example of the procedure: * = Lu] +N [u] +f [1] 


The operator Lu] is linear, and N [u Jis a nonlinear one. Both in the linear and in the nonlinear 


operator appear the field u(t, x) and its space derivatives, multiplied by suitable coefficients that can 
be functions of time and space coordinates. We assume that in the nonlinear operator there are no 
linear terms (all linear terms appear in Lu] by construction). The term /f is a source or external 
forcing function f (t, x); 

Furthermore, we suppose that the zero solution is our reference solution of equation [1], although in 
the original problem we may be faced with a steady state field with spatial variations. If this is the 
case, the equations are manipulated until an equation like [1] is obtained (this can be always done). 


Then we can expand the nonlinear operator in a Taylor series, using the Fréchet derivatives 
introduced in functional analysis (Griffel, 2002): 


Nu] =5- D*[osw,u]+ =-D*[O:u. uu]. [2] 


The operator D* [o; Uu, u|is bilinear, D* [o; Uu,u, ul is trilinear, and so on. 


We assume that the linear operator has simple eigenvalues A , with eigenfunctions @g, such that the 
set of eigenfunctions is complete in the space of functions that has the considered field, at each 


instant, as one of its elements. 
Now, let us consider the adjoint operator L’ [u] corresponding to Lu], its eigenvalues 1, and its 
eigenfunctions Q p (z): L lo’, |= A, ‘Dp [3] 


The set of eigenfunctions of the adjoint operator is also a complete set in the function space to which 


the field u(t ,x) belongs. Each eigenfuction 2, is orthogonal to the eigenfunctions g, of Liu|that 


corresponds to the other eigenvalues: 
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It is always possible to put: (9; ; ?,) =1 Next, we define the mode amplitude A z (t) by: 
AylO)=(9" vou) = [ nl) ule. ¥)- wl) (5] 


Here w(X) is a positive weight that may be identical to 1. Then, with the modal amplitude A " (t) 


fo) 


given by equation [5]: u(t, x) = > A, (t)- 2, (x) [6] 


p=0 


Making the internal product of both members of [1] with ?, we obtain: 


(0) = (05 tu) + (95. Mlu) + (95.2) im 
300 Apc dA, (t) 
But (0; = = rage ) = a; [8] 
(9;-Lul) = (cp; }u) =2, -4,(0) [9] 
(¢}.Mlu)) =5-(9,.P *Io; Uu, u) + =-(9,.D'[Osu.u.u) + [10] 


Now we substitute the expansion [6] in equation [10], and after several calculations it follows: 


(7, ul) = yo P.qr 2A, “A, = YO, 0. iD “A, uA, Tee [11] 


qg.r=l q.r,s=l 


Here, by definition: 
0,4, =(9;.D"|0:9,.9,]) 2a 9, are =(P2.D'[0:9,.9,.9,]) 11201 


7] 


Some of the symbols 0 eres 


one . as well as the eigenvalues A, may be functions of time. 


Putting together equations [7] to [12] and defining the projection of the source term onto @"» (x) 
by f, (t) = (9; my i ) , we derive a numerable infinite set of ordinary differential equations in the mode 


amplitudes: 


2 =A,-A .+ ye ar 4, A, + yo vars 4, A, A, t+ f(t) p=2... [13] 


Ee q.r=l g.r,s=l 
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In practice, a truncation must be done, so that a finite set of nonlinear ordinary differential equations 
remains. In the linear approximation all the modal amplitudes appear uncoupled: 


dA,(t) ere 

Ta =A, ‘A, (t)+ ‘A “ (t) As the superposition principle does not apply in the case of 
nonlinearly coupled mode amplitudes, one may ask about the utility of nonlinear modal analysis. The 
nonlinear interaction consists in the exchange of energy between different modes. Thus, even if the 
initial field may be represented by only one or a few modes, the nonlinear processes tend to activate 
other modes as time progresses. In a Fourier setting this would be a generation of harmonics or sub- 
harmonics. Besides the possibility of creation of spatial patterns of smaller or larger wavelengths, 
there are other very important effects: self-damping of a growing mode and inter-mode suppression. 
The self-damping may appear in a mode that is unstable according to the linear theory 


(because A, > 0). If this mode is the only one that is perturbed from its steady state its amplitude first 


begins to grow exponentially: A, (t) =A, (0) : exp(Z, : t) When @ <0, the quadratic term tends 


P>P>P 


to slow down the growth of the amplitude (self-damping). 


If the initial perturbation affects several different modes and if initially many of these grow according 
with the linear approximation, the growth of one mode could tend to damp the growth of the other 
modes, either because it is the fastest growing mode or because its initial amplitude is much larger 
than the initial amplitudes of the other modes. This occurs when there is competition between modes, 
for example a competition for a given amount of energy. Both self-damping and inter-mode 
suppression allow a partition of the space of mode amplitudes into basins of attraction of equilibrium 
points, leading to bi-stability, run-away, and other interesting threshold dynamic behaviors. 


When there are dominant eigenvalues, the evolution equations cast according to the above mentioned 
summarized restricted procedure of nonlinear analysis are suited to study mode slaving in Haken’s 
sense and to apply Eckhaus’s method to determine sufficient conditions for instability. 


Applying suitable asymptotic methods, analytical formulae can be derived to characterize the dynamic 
behavior of a small number of dominant modes, and the understanding thus obtained serves as a guide 
to the design of in-depth digital simulations, using modal equations or directly by discretization of the 
original nonlinear field equations. 


Nonlinear modal analysis of the cardiac bidomain equations. 


1. Bidomain equations: 
To describe the polarization of myocardium membranes due to electric current field produced by an 
external electrode, it is possible to recast the equations of the mathematical model of bidomain as 
follows: 


ae (cu EV (.7)+ Soo Var i))=velEevv, nv (5-7) [14 a] 


Ww =FV,.W) [14] 
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Here f represents time instant, 7 represents the position vector in three-dimensional Euclidean space 
and vectors are represented with an arrow over the corresponding symbol. The wavy bar over a 
symbol represents a tensor (dyadic) and V represents the gradient operator. The gross dot ¢ 
represents a scalar product of two vectors of a tensor (dyadic) with a vector or of two tensors. As 
usual V o ( ) represents the divergence operator. The small dot - represents the product with a 


number. 


The electric current density J (t,7 ) in the volume conductor of the tissues, that appears in the last 


term of the second member of equation in [14 a], is assumed known. 


The continuum of excitable membranes is described by the state variables (Vv, (t, F ), Wt, r )) where 


Vy (t, r ) is the field of membrane voltages and Wt, 7) is the vector field whose components are 


channel activation and recovery variables. The time evolution of these channel variables is given by 
equation [14 b]. The ionic current density through the excitable membranes is represented 


by Ji, Vv, W), C,, is membrane capacity per unit area and y,, is the area of membranes per unit 


volume of the bidomain. The tensors that appear in the first and in the second terms of the right 
member of equation [14 a] are given by the formulae: 


C=f,-G,eG'ef,-G, [l4c] b=>. f,.G,-f,-G, eG [14d] 


The scalars f; and f, give the volume fractions of the intracellular and extracellular continuum, 


respectively, while the symmetric and positive definite tensors G, and G, represents the electric 
conductivity of the intracellular and extracellular continuum, respectively. Here, by definition, 


G=f,-G,+f,-G, is a positive definite and symmetric tensor that characterizes the electric 


conductivity of the heart tissue, and Gis its inverse tensor. From formulae [14 c] and [14 d] it 


follows that both tensors, C and D are symmetric and positive definite also. 


The coupling between the electric current field in the tissues volume conductor and the polarization in 
the continuum of excitable membranes is given by the forcing term V e (5 eJ ) in the right-hand 


member of equation [14 a]. This forcing term is called a generalized activating function (Tung, 2009). 
For the quasi-stationary fields of interest in cardiac pacing the field of electric current density is 
divergence free: V J =0 Then, if the intracellular and extracellular continua are homogeneous and 
isotropic, the generalized activating function vanishes identically. In this last case, as will be seen in 
section (B), the coupling between the electric current field due to the electrode and the excitable 
membranes appears through the boundary conditions at the interface that separates the myocardium 
from the other tissues. If the intracellular and extracellular continua are isotropic and heterogeneous, 


the tensor D is a scalar function of position D(7).Then the activating function reduces to 


VD( r )e J (t, ii ) and in a non-zero current field it vanishes only if the current density vector is 


orthogonal to the gradient vector. If the intracellular and extracellular continua are homogeneous and 
anisotropic, the coupling term vanishes identically only if the anisotropy is equal. But if it is unequal, 
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as it is in the ventricular myocardium, the activating function doesn’t vanish, and the coupling that 
results produce the so-called virtual anode and virtual cathode effects (Tung, 2009). 


As the external current fields are quasi-steady, the current density can be expressed as the product of a 


spatial field F( 7 ) by the electric current i(t) injected by the electrode: 
I(t,7)=i(t)-F(7) [15] 
Then the generalized activating function is factorized thus: 
V 0(De0J)=i(t)-V0(De F(7))=i(t)-K(7) [16] 
The following scalar field will be called form factor of the generalized activating function: 
K(7)=Ve(D«F(7)) [17] 


2. Boundary conditions, a constructive definition of regions of influence and 
identification of a suitable linear operator 

Let ai(7 ) be the field of unit normal vectors at each point of the interface between the myocardium 

and the outside tissues. Assuming that in the bidomain’s boundary the electric current goes to or 

comes from the surrounding tissues only through the extracellular domain, in the surface of the 

cardiac muscle the current density in the intracellular continuum must vanish. If the electric current 

density in the intracellular continuum of the myocardium is supposed to vanish at the interface, it is 


possible to show that ne (C eVVi )= ne (, G, ® (G" ° i) [18] 


Equation [18] couples the polarization of the continuum of excitable membranes with the electric 
current field produced by the pacing electrode at the interface between the myocardium and the extra- 
myocardial tissues. It doesn’t vanish identically, even if the intracellular and the extracellular continua 
are assumed to be homogeneous and isotropic, so boundary condition [18] is the only always present 
coupling mechanism between the electrode and the excitable membranes. For a quasi-steady field, 


equation [18] reduces to: ne (c eVVi )= n o(, G, ° (G" e F(7))). i(t)= i(t)- S(7) [19] 


The scalar field S (7), defined in the interface between the myocardium and the extra-myocardium 


tissues, is a surface form factor: 

5(7) =relf, G, (G+ F())) [20] 
When the active electrode is located near enough to the myocardium, the polarization of the excitable 
membranes is determined mainly through the spatial variation of the applied external electric current 
field (the field F ( e ) in equation [15]). This variation acts on the excitable membranes through two 
mechanisms related with the applied field. One is the form factor K( r ) (given by equation [17]) of 


the generalized activating function. Let us call this the volume form factor. The other is the surface 
form factor S (7 ) that appears in the right-hand member of the boundary condition [19]. 
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Now, a region of influence of the active electrode in the myocardium must be defined, to use it as 
the bounded spatial domain B needed to develop a nonlinear modal analysis. For this purpose, let us 
consider again the surface and the volume form factors. The surface form factor is a scalar function of 
the position in the boundary surface of the bidomain. The volume form factor is a scalar function of 
the position in the interior of the bidomain. Now, consider the family of curves (in the boundary of the 
bidomain) or surfaces (in the interior of the bidomain) in which the surface form factor and the 
volume form factor, respectively, are constant. 


Next, consider the lines orthogonal to the curves or surfaces, and take the restriction of the 
corresponding surface or volume form factor to one of these lines. We thus obtain the values of the 


form factor K(F) or S(F) as a function K,(s) or S,(Ss) of the arc length s corresponding to the line 
I’ that we are considering. After that we follow the line, leaving the region that is significantly 
polarized. We will arrive at a final inflection point s,; of the functions K,(s) orS,(s), after which 
the corresponding form factor will tend monotonically (and usually quickly) to zero. Then, we 


calculate the following estimate of the distance (along the orthogonal lineI’) in which the absolute 
value of the form factor suffers a significant decrement towards zero: 


pox eo [21 al as= Bel [21 b] 
i K;(s;) GS) 


Next, we identify the point of I” whose arc length is given bys, + As, and we repeat the operation 


for each possible orthogonal line. Finally, relative to the volume form factor, we can define the 
boundary of the region of influence as the set formed by these points, and the region of influence 
B as the bounded set of points inside this boundary. Thus, we have the bounded domain needed to 
develop a restricted nonlinear modal analysis of bidomain models. 


The boundary OB of the region of influence can be split in two surfaces: OB, and OB, 
While OB,is located near the electrode, at the interface between the myocardium and the other 
tissues,OB, is located inside the myocardium, in the region of transition between the polarized and 


the non-polarized membranes. 


The following non-homogenous Neumann boundary condition is applied ind B, : 
iie(Cevv,, )=i(t)-S(7) [22] 


If V, is the rest value of membrane voltage and W,, are the rest values of the channel variables, a 


homogenous Dirichlet boundary condition is applied ind B, : 


Vy =Vp [23 al W,, =W, [23 b] 


Now, it is necessary to identify the linear operator whose eigenfunctions and eigenvalues will be used 
in a nonlinear modal analysis. 
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Equation [14 a] can be rewritten thus: 
O = e 
ty 5 Vult?)=-Ry idee (V,,.W)+Vve(Z,? eVV,, Js i(t)K,(7 ) [24] 


Here R,, is a suitable unit area resistance of the membrane, T,, = Ry -Cy is a membrane time 


ae eae : de Ses 
constant, 2,,, =—“-C is a tensor space constant and K,(7 )= v«[ de () is a 
Xm Am 


modified volume form factor. The equations [14 b] of evolution of the channel variables remain 


unchanged: SW = Flv, Ww) 


The ionic current is developed in a Taylor series relative to the rest state: 


ov, V, W a (Vi. -Vp)+ +e (V,.W, Jo (Ww -W,)+ higher order terms 


Jy, Vy .W) = 


Then, the following linear and self-adjoint operator is introduced: 
LV yp ]=—-Ry 2 (Vp We) Vue —Ve + V 8 (Ay? © VV | [25] 


It operates on real regular functions defined in the influence region B. These functions verify 


homogenous Dirichlet boundary conditions ind B, and homogenous Neumann boundary conditions 
inOB,. Membrane voltage and channel variables are developed in series of the normalized 


eigenfunctions of the operator: 
1 Ve DA- 066 W,, =, +E, (7)-B,(0) 
q 


To find the equations for the mode amplitudes, following a procedure like the one summarized in Part 
2, both J, (V,, Ww) and F V,, Ww) are approximated by Taylor polynomials in terms of the 


deviations of the state variables from their rest values, beginning with a suitable excitable membrane 
model. 


3. Nonlinear modal equations for the cardiac bidomain with a simplified mathematical 

model of membrane’s ionic current 
If the activation variables are relaxed to equilibrium with membrane voltage, and a single recovery 
variable is considered, besides the well-studied two variables membrane model of Fitzhug-Nagumo, 
there are simple mathematical models specifically constructed to describe threshold effects in cardiac 
tissue (Aliev and Panfilov, 1996) that could be a suitable starting point. However, to simplify the 
analysis, the recovery variable will be considered frozen. The only membrane state variable that 


remains in this case is V,,. 
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To capture thresholds effects, order three in b = V,, — Vz is enough for a polynomial approximation 
to... (Vi, )= nee (v). So, let us consider the following formula for the ionic current density, with b 
and c always positive: Ry ‘Tog (V ie —v+b-v*-c-v’ [26] 


The linear Sturm-Liouville operator [25] can be written Lv] =Ve (7 eV v)+ D 


Its orthogonal eigenfunctions defined in the region of influence of the electrode and satisfying mixed 
Dirichlet and Neumann homogeneous boundary conditions there, can be determined solving: 


-Veli, eVo,)= U, -®, pal 2c [27] 
Following the procedure outlined in Part 2, with A» = (i + = and employing normalized 


eigenfunctions, we obtain for the mode amplitudes: 


dA ee ce) 
to =A,-A,+b- Lie -A,-A,—c: ae -A,-A,+A, +..+ f(t) 
dt q.r=l q.r,s=l 
pel 2sc [28] 
Here, by definition: 
Tar = [Pp Pq°P AV [29a Lyars =] Pp ‘Py Pr QsdV [29] 
B B 


Basel - 
If S, (7) =—“.$§ (7) the forcing terms are given by the formulae: 
Xm 


OB, 


1) Je, (0) )-dV+ [e,(7) le) p=1,2,3,... [30] 


The volume forcing term fav oh | 9, (F K(P)av) gives the projection of the volume 


form factor (of the generalized activating function) onto the eigenfunction 9,(F), and allows the 


study of virtual anode and virtual cathode effects on threshold dynamics. It vanishes when the 
activating function vanishes. 


The surface forcing term f, (t)= i(t): | 2, (7). sea gives the projection of the surface 
form factor onto the eigenfunction @,, (7 ), and is always present and operating to polarize the tissue 
during an applied current pulse, even in a homogeneous and isotropic syncytium. 

4. Alibi, eigenfunctions and eigenvalues 


Now, let us focus on the eigenvalues and their corresponding eigenfunctions. For common types of 
heterogeneity and anisotropy, like the ones found in ventricular myocardium, a change of spatial 
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variables from Euclidean to curvilinear coordinates may be found, allowing the reduction of a 
complex three dimensional membrane polarization problem to the polarization of membranes in a 
homogeneous and isotropic case, if an auxiliary “virtual space” (alibi) with a smooth one to one 
correspondence with physical space, is introduced (Suarez-Antola, 1994 and 2008). It is possible to 
apply this idea of an auxiliary space (alibi), to the description of the polarizing effect of an active 
cardiac pacemaker electrode. In the auxiliary space the eigenfunction-eigenvalue problem easier to 
solve by analytical methods. Going back to the original coordinates, the results are interpreted in 
physical space, and the modifications that the polarization of the excitable membranes suffer due to 
the anisotropy and heterogeneity of the tissues can be assessed. 


We assume that a small active electrode is in the ventricular cavity, touching the endocardium. The 
endocardium surface may be represented, approximately enough, as a plane. The myocardial fibres 
run, almost parallel, in layers. However, the direction of the fibres in each layer varies smoothly from 
the endocardium to the epicardium. 

We introduce first a trihedral T, of unit vectors, ie as k, with the origin and the unit vectors i. I 
fixed in the endocardium plane, and k, orthogonal to this plane, pointing towards the epicardium. 
Taking as new origin each point of coordinate z of the half straight line defined by the unit vector k, 
and the origin of T, , we introduce a new orthogonal trihedral T, of unit vectors é, = é;,€,,€; = ko, 
where é y gives the direction of the fibres in the layer located at a depth z from the endocardium. Then 


any point in the myocardium may be given by the coordinate z and its three coordinates 


X,,X,X3 relative toT.. 

The following connection exists between the coordinates x, y relative to T, and x,,x, relative to T, : 
Xx, =Xx- sin[a(z)|+ y: cos[a(z) [31a] 
X, =X: cos[a(z)]+ y: sin[a(z)| [31b] 

Here a(z) is the angle between the unit vector é , (that gives the direction of the fibers in the layer) 


and the unit vector - : 


cs Ro 
Considering that the tensor space constant is given by the formula ye = —“ .C, equation [27] may 
M 


| ae 
be written thus: Ve -Ce ve, | = Ul, nG [32] 
Xm 


Due to the symmetry of the intra G, and extracellular G, conductivity tensors (Antoni, 1998; Janks 
and Roth, 2009; Henriquez and W. Ying, 2009; Tung, 2009), it can be assumed that the tensor 


C= f, -G, eGie f, -G, is diagonal relative to T, in the corresponding layer of myocardial fibres. 
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With enough approximation, its diagonal components C,,,C,,,C,, can be represented as functions of 


only one coordinate, that is C,, (x, ) i oe (a ) and C,, (x, ) : 


Therefore, if @, (3 X>, x) is the representation of the eigenfunction in local coordinates, equation 


3, O(R O 2 
32] reduces to: apne eae ee ae 2 
[32] reduces to 2. ae 0) ax. °,| Hy “Pp [33] 


Nevertheless, we are still in physical space, naming the same points in another form (alias, using 
other coordinates). 


Now, let us introduce a virtual space (alibi) of coordinates €,,¢,,¢,, filled by a homogeneous and 


R 


isotropic electric syncytium. The scalar space constant of this auxiliary syncytium is 2,, = ,|—“-C, 
M 


, with a scalar C, in place of the tensorC . The connection between corresponding points of each 


ae * | C 
space is given by: Gee J C6) ‘ds kk =1,2,3 [34] 


The connection between the local representation of the eigenfunction in physical space 


2, (x,,.x,,x, )and the eigenfunction in the corresponding point of virtual space virtual space 


y ,(&,.£).€,)is defined by: V (G1...) = @, (a (E), 2(E ), 23(&)) [35] 


Furthermore, we introduce the functions: Cy (é, ) =C,, (x, (é, )) k =1,2,3 [36] 


R 
From — - C, = ys , [34], [35] and [36], the left-hand member of [33] may be shown to be equal to: 


M 
3 @ 3 Fa) 
24 ae 7D) te) (22) 
k=l k Cy dg, Sk 


ae L 
So [33] is equivalent to: Ag ey, me C = vt —wW, [37] 


Here: Ly x ® [38 a] Lek © Rt [38 b] 
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If the spatial scales of the components of the gradient of the eigenfunction in virtual space L,, ,, are at 


least an order of magnitude less than the spatial scales of the components of the tensor C in virtual 


space L,,, , then the eigenfunction in the auxiliary space y, (é, Gis é,) will verify, approximately 


OW 2 
the equation: Au oF 2 zw W, [39] 


The approximation will be better the smaller are the quotients of the spatial scales Ey, i ONL, 4 


Is seems reasonable to assume that the region of influence of a symmetric electrode over the 
membranes of the myocardium cells in the alibi space may be represented by the half ellipsoid of 


revolution, with axes R, and c, shown in Figure 1. 


Figure 1 Region of influence in alibi space 


2 2 2 
+ 
The equation of this half ellipsoid is: & z $2 + $3 zy =l fy 20 [40] 
0 Co 


In this formula R, is a certain radius of influence of the electric current field over the surface of the 


[41] 


(alibi) electric syncytium, and: 
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(a, = 5,783) is a measure of the depth of penetration of the influence of the active electrode in the 


(alibi) electric syncytium’. 


The eigenfunctions y, are Lamé functions obtained after writing equation [39] in ellipsoidal 
coordinates (Polyanin and Manshirov, 2007), using the resulting symmetry and applying 
homogeneous Dirichlet boundary conditions on the curved surface of the half ellipsoid. From 
equations [34] and [35] the eigenfunctions @,,(x,,X),X3)in physical space can be derived. The 
fundamental normalized eigenfunction @, doesn’t change its sign in the region of influence. 
However, when the index p increases, the number of sign changes in the normalized eigenfunctions 
g,in the region of influence increases. This result follows from the general properties of the 
eigenfunctions in a Sturm-Liouville problem (Stakgold and Holst, 2011). Then, from [29 a] and [29 b] 
we can see that the coefficients J pqr and I pars decrease when any index increases enough, some 
being positive while others being negative, and others equal to zero, depending of the behavior of the 
product of eigenfunctions in each integrand. 


The eigenvalues can be estimated by the following approximate formulae (for a heuristic derivation, 
see the references in foot page note 1): 


2 
5-7 A 
el =| i+ b+a,"){ 4 p=1,2,3.... [42] 


0 


In [42] the numbers @ pare zeroes of the Bessel functions of the first kind and order zero, so 
that a, increases quickly with p: @ - =5.7830, @ a = 30.4715, a = = 74.8865, ... From [42] 
we see that the negative parameters A, = (1 +u 2) that appear in the coupled evolution equations 


A 
for mode amplitudes (equations [28]) take values widely separated if —“ is of the order of 
0 


magnitude of | or higher, and can be ordered thus: O> A, > A, > A, >... [43] 


5. A modified Eckhaus method 


Next, let us estimate the size of a perturbation in membrane voltage which causes instability. To apply 
Eckhaus asymptotic method, two restrictions must be satisfied (Eckhaus, 1965; Denn, 1975): 


(1) The shape of the projection A,(0)= | v(0, r )- Q, (7): dV of the initial disturbance onto the 
B 
first eigenfunction must be dominant over the other projections. 


“ The complete derivation of the region of influence over a homogeneous and isotropic electric 
syncytium, as well as approximate formulae for the eigenvalues and eigenfunctions, can be found in 
Suarez-Antola, 1994 and 1999. Some consequences of these and other related results, from the point 
of view of for cardiac pacing, are briefly discussed in Suarez-Antola, 2006. 
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(2) The dynamics can be studied near a point of marginal stability where A, TO, 


These restrictions don’t apply to the case of the threshold dynamics of excitable biological tissues 
stimulated by external electrodes. Nevertheless, it is enough to have separated enough eigenvalues 
and to have the shape of the projections of the initial disturbance onto the first eigenfunctions 
dominant over the other projections, to be able to study threshold dynamics with a few low order 
mode amplitudes (Suarez-Antola, 1994; Korenko, 1997; Suarez-Antola, 2012). 


6. The global features of threshold dynamics in the space of mode amplitudes: 
decaying sets, amplifying sets, threshold manifolds and reachable sets 


To consider virtual anode effects during cathode make excitation, as well as to study cathode break, 
anode make and anode break excitation, considering virtual electrodes effects, it is necessary to work 
with higher order mode amplitudes and at least a recovery variable in the mathematical model of the 
excitable membrane. For symmetry reasons, mode amplitudes of even order are often not perturbed 
from the rest state. 


So, it may be enough to work with the first two or three odd mode amplitudes as in the case of the 
study of threshold dynamics of a nerve fibre stimulated by external electrodes (Korenko, 1997; 
Suarez-Antola, 2012). However, a global qualitative pattern emerges, common to all cases of 
nonlinear modal analysis of threshold dynamics. Figure 2 sketches the main aspects. 


Amplifying Set 


Reachable 
Sets from 
the origin 


Pseudo-Threshold manifold 


Figure 2 Qualitative patterns of threshold dynamics in the space of mode amplitudes 


For each electrode-tissues-target element system, the space of mode amplitudes can be approximately 
divided in two subsets: the decaying set and the amplifying set. The rest state, represented by the 
origin of the space of mode amplitudes, belongs to the decaying set, as well as those points such that 
if the system is there at the end of the stimulating pulse, the system will return to rest without showing 
an action potential. If the pulse leaves the system in the amplifying set, an action potential will 
emerge, and the excitation succeeds. The states located at boundary between the decaying and the 
amplifying set can be considered as threshold states: for each electrode-tissues-target element system 
there is a whole family of these states. When the recovery variables are considered, the boundary 
between the decaying and the amplifying set is fuzzy, so in this case there is a pseudo-threshold 
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manifold, as suggested in Figure 2. Besides, this figure shows the reachable sets whose states can be 
attained by applied current pulses beginning from the rest states of the excitable membranes. 


Results for cathode make excitation described with one mode amplitude. 


For a perturbation due to a cathode, such that the first mode amplitude can be taken as dominant over 
the other mode amplitudes, it is possible to study threshold in a first approximation to work with the 
equation of the first mode amplitude uncoupled: 


fy he (L 4g VA, + Dy Ae = egy Ay +f (¢) [44] 
f(t)=i(t)-| [o,)- Ko(#)-av + [ o(7)-5,(7)-as [45] 

B eB, 
> (5° | 1 147-7830.{ 4u [46] 

ae R, 


For a convex cathode, the radius R, of the influence region in alibi space can be estimated, as a 


function of electrode radius 7, and the thickness e of unexcitable tissues interposed between the 


electrode and the excitable myocardium, by the formula 
R, =7-(n +e) [47] 


In this formula the proportionality factor vy considers such things as monopolar and bipolar electrodes 


and local conductivity variations associated with inflammatory responses to foreign body (Suarez- 
Antola, 1994 and 2006). 


1. Thresholds: liminal regions and families of threshold states 


Let us suppose that a current pulse injected through the cathode and ends whent =0. The pulse 
carries the mode amplitude from its rest value A, =O to a certain value A, (0)> 0. After this, the 


mode amplitude evolves according to the free dynamic equation: 


dA 
tM a FH AL + By Ay Clu Ay [48] 


The equilibrium amplitudes of [48] are shown in the bifurcation diagram of Figure 3 as functions of a 


parameter P o u, 


2 
When the discriminant , -( pli ! (1 fb ) is negative, the only equilibrium amplitude 
2¢liin I 


1111 
corresponds to the rest state A, = 0 of the excitable membranes: in this condition the syncytium is not 


excitable. 
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When the discriminant is positive, there are three equilibrium amplitudes: 


b-I 
A,=0 [49a] A, = AY al A [49 b] A= 
2cLiy 


bly 
2cliiny 


+VA [49 c] 


The lower branch in Figure 3 gives A , as a function of P o 44, , and its points correspond to unstable 


(threshold) equilibrium states. The upper branch gives A , , as a function of the parameter P oc 4, . 


le 


Figure 3 A global bifurcation corresponding to one mode threshold dynamics. 


2 2 
b Ty 2 ir 


4-¢? 


The discriminant A =0 for acritical value jz, , that verifies1+ 7 .= 


When //, < //,,,the decaying set is given by the inequality A, < A the amplifying set is given by 


lu? 


A, > A, and the threshold manifold reduces to the single point A 


Lut 


5-7 to 
From equation [46] and [47]: Mt, = .|| —— |-| 14+ 7,7830-| —_—. [50] 
8 y° (7, as e) 


As consequence, the threshold amplitude A besides being a function of the parameters of the 


lu? 
excitable membranes, is a function of electrode-myocardium distance % +e. As this distance may be 


different for different pacing electrodes implanted in the same heart and may vary for the same 
electrode during the so called “physiological threshold increase”, the threshold of cardiac pacing is a 
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systemic property. As remarked in section (F) of Part (3), there isn’t a single threshold, but a family of 
threshold states. 


The restriction “/, < 44, ,must be verified to be able to elicit an action potential. From equations [46] 


and [41] this restriction can be expressed as a restriction imposed to the radius and depth of the region 
of influence of the electrode in alibi space: 


Ryo R= 2D Aa [51 al 
8 2 
5 : nr : 1. = 
R 
Co > Coe = = [51 b] 


R 2 
|s30 + [Fe] 
Au 


If the inequalities [51] are not verified, the excitation will fail. Going back to physical space, 
nonlinear modal analysis predicts the existence, shape, and size of liminal regions in the sense of 
Lindemans (Lindemans, 1977). 

The region of influence of the electrode must include a liminal region to be able to elicit a propagated 
action potential. 


2. Time constants and rheobases for current pulses 


Taking the linear approximation to equation [44] in the decaying set, where A,(t)< A,,,, and re- 


ordering it a little, it follows: es x Es A, + 2, g, u(t) [52] 
dt is Tai 
Here it appears a time constant defined by: t= A fu ; [53] 
+h, 


The combined projection of the form factors onto the eigenfunction of the assumed dominant mode is 


given by: C= [o(F)- K,(7)-av + [a(F)-s.(7)-4s [54] 
B 


oB, 


Integrating [52] between tf = O and ft =d (beingd the duration of a cathode pulse), and assuming that 


the excitable membranes were initially at their rest state (A,(0) =()), at the end of the current pulse 


d 
the dominant mode amplitude verifies: A, (d ) = | G(d - t) . i(t)- dt [55 a] 
0 


G(t)= = f, on - a [55 b] 
™ f, 
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A current pulse of duration d_ will be just threshold if A (d ) = A,,,. Applying this last condition, 


jointly with [54] and [55], to a rectangular current pulse the Lapicque-Hill-Blair relation between 


pulse amplitude I, (d ) and pulse duration at threshold is obtained, as a byproduct of nonlinear modal 


analysis and the linearization hypothesis in the decaying set: 


I 
I, (@) = ——2_ [56] 
1— exp} — — 
f, 
l+y,)-A 
Current rheobase is given by: Ie = Mae ead [57] 
§1 
Limit threshold charge Q, , is given by: QO, =lim,,, (d 7 ee (d)) =Il, acts [58] 


From equations [50] and [54] together with [49 b] that relates 4, and the threshold amplitude A 


ee 
the time constant (proportional to the chronaxie), the rheobase and the limit threshold charge are 
obtained as well defined analytical functions of geometric and bioelectric parameters of the electrode- 
myocardium system. The predicted increase of the time constant (chronaxie) and rheobase when the 
electrode radius and the thickness of interposed non excitable tissue increases, agrees with known 
experimental results (Lindemans, 1977; Irnich, 1975, 1980 and 1983; Laplantif, 1983). 


a. Time constants and rheobases for voltage pulses * 


Modern permanent pulse generators use constant voltage pulses to pace the heart. While the voltage 
remains at the programmed value, the electric current suffers variations that reflect the changes in the 


load impedance seen by the pulse generator. This impedance is the sum of an ohmic resistance R,, and 


the impedance produced by the polarization of the electrochemical double layer located at the 
interface between the electrodes and the tissues. The voltage drop between the anode and the cathode 
U (t)can be written, in a linear approximation, as an ohmic voltage drop in the cables and the tissues 


is 
plus the voltage drop in the double layer: U (t) =R,1 (t) + | K (t = u)I (u).du [59] 
0 


Here K (t)is positive and | K (u).du =R,, being R, the D.C. resistance of the interface. As U (t) is 
0 


given by the voltage pulse generator, it is convenient to invert [59], thus obtaining: 


=z ve) [uc “ula (60) 


i, R 
The kernel L(t) is also positive and | L(u).du =—__., 
‘ R,+R, 


For a voltage pulse of constant amplitude U,, , equation [60] may be rewritten thus: 


5 This section is taken from Suarez-Antola, 2006. 
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I(t)=——*—- - (1+ lt [61] 
(= gag (+96) 
lim t U, \d 
The fu8nction y(t) is positive, monotonically decreasing from ft, ., = aol e ral ) towards 


Up 


zero whent > ©. 


t 

It can be given as a function of } L(u)du. From [61] it follows that / (t) decreases from the beginning 
0 

to the end of the voltage pulse. 

According to the analysis done in the previous section, a voltage pulse of constant amplitude will be a 


threshold pulse U, , if the corresponding current pulse J, t verifies: 
d 
A, d AiG d—t-I,t -dt=A, [62] 
0 


Substituting the threshold current pulse given by Equation [61], written for the amplitude i of a 


threshold voltage pulse, in [62], we obtain: 
(R, +R,). Oro 


Uz o(@)=-—_+ [63] 
[aw at+ | G(d —t)-g(t)-dt 
0 0 
ait d 
IfG(t) =e '’ , itis possible to show that: lim,,_,,, [ca - t)-o(t)dt =0 
0 
But } G(r). dt =t,, so considering that/, = ay, , the rheobase of voltage is: 
5 s 
U, =lim,,, Uy o(d)=(R, + Ry) Te [64] 


Voltage rheobase and current rheobase are connected by the total resistance (D.C. impedance) seen by 
the pulse generator. 


As the parameter f,is the time constant of the strength-duration curve obtained with controlled 
current, from now on it will be represented by ¢ Sue 


The strength-duration formula for rectangular voltage pulses, given by equation [63], can be rewritten 
thus: 


be Un 


U,, d = [65] 


d d 
eae eee 


The strength-duration formula for rectangular current pulses obtained from Equation [62] may be 
rewritten thus: 


Iho d a [66] 


For short pulses they both behave proportional tod but for longer pulses the threshold amplitude of 
rectangular voltage pulses approaches to its rheobase faster than the threshold amplitude of 
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rectangular current pulses of the same durations. All this is in accord with experimental results 
Urnich, 1975 and 1980; Furman, Hayes, and Holmes, 1993). 


lim, (dU 76 (d ) 


From Equation [65], considering that by definition?, |, = 7 
R 


, it follows: 


R 


t = ts 7 = Pp .t 
ae EEO) Rea, 


Thus, the time constant for excitation with voltage pulses is always less than the time constant for 
excitation by current pulses, and the difference between them increases with the increase in the D.C. 
impedance of the electrochemical double layer relative to the ohmic impedance of the lead cables and 
the tissues. All this is in accord with known experimental results (Irnich, 1975, 1980 and 1983; 
Furman, Hayes, and Holmes, 1993). Figure 4 summarizes some results. 


[67] 


ms Sl 


- 


s,U 


Figure 4 Relation between the time constant and the electrode radius for controlled current and 
controlled voltage rectangular pulses. 


Conclusions: 


Several experimental results are explained analytically: the existence of liminal regions and the 
behaviour of the parameters of strength-duration curves as functions of electrode-tissue distances, 
electrode size and cardiac bidomain parameters for cathode make excitation. 


Some new facts are predicted analytically: the size and shape of the liminal regions considering 
tissue’s anisotropy, the family of electrode-myocardium patterns of threshold states, and the 
behaviour of the parameters of the strength-duration curves as given by the formulae obtained 
applying nonlinear modal analysis for cathode make excitation. 


In principle, nonlinear modal analysis could be applied to both, cathode make, and brake and anode 
make and brake excitation of myocardium, considering virtual cathode and virtual anode effects by 
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means of equation [30] that gives the projection of the generalized activating function onto the mode 
eigenfunctions. 


The main shortcoming of the proposed general procedure is that it is unable to cope with the growth 
of the action potential, being restricted to just-threshold dynamics. However, from one side this may 
be enough for the discussion of subjects of interest from the viewpoint of clinical cardiac 
electrophysiology, while from the other side the obtained results could be used as guidelines for doing 
both, new experiments, and new in-depth studies of digital simulation of myocardium threshold 
dynamics. 


The method outlined here can be applied, in principle, to the excitation of nerve, skeletal muscle, 
myocardium and smooth muscle, not only when the active electrode is convex but also when it is 
concave. It is possible to determine theoretically the form factor produced by a concave electrode in a 
functional syncytium (Suarez-Antola, 1994 and 1999). Then, some of the observed differences in 
behaviour between strength-duration curves produced by convex and concave electrodes used to 
stimulate the same myocardium (Suarez-Antola, Griego and Fiandra, 1994) can be at least 
qualitatively explained. 
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